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FOREWORD 


This  final  report  wa3  prepared  by  the  Structural  Mechanics  Section  of  the 
Grumman  Aerospace  Corporation,  Befchpage,  New  fork,  for  the  Vehicle  Dynamics  and 
Structures  Division,  Air  Force  Flight  Dynamics  Laboratory,  Wright-Pat terson  Air 
Force  Base,  Ohio.  The  work  ’as  performed  under  Contract  No.  F33615-72-C-1101, 
which  was  Initiated  under  Project  No.  1370,  "Dynamic  Problems  in  Flight  Vehicles", 
Task  No.  01,  "Aeroelastic  Problems".  Initially  Mr.  R.  F.  Taylor  (ITS ) and 
Dr.  V.  B.  Venkayya  (FBR)  were  the  Project  Monitors  of  this  contract,  after  which 
Capt.  S.  M.  Batill  (FIS)  assumed  this  position. 

The  report  consists  of  two  volumes.  Volume  I, entitled  "Theory  and  Applica- 
tion", describes  the  analysis  and  redesign  procedures  provided  by  a computer  pro- 
gram system  for  minimum-weight  design  of  cantilever  or  free-free  lifting-surface 
structures  subject  to  combined  strength  and  flutter-speed  requirements.  Detailed 
instructions  required  to  use  this  Flutter  And  STrength  Optimization  Program 
(FASTOP)  are  provided  in  Volume  II,  entitled  "Program  User's  Manual".  The 
report,  which  covers  work  conducted  between  1?  March  1972  and  31  December  1975, 
was  submitted  to  the  Air  Force  in  December  1975. 

Dr.  W.  Lansing  was  the  Program  Manager  and  Mr.  K.  Wilkinson  was  the 
Project  Engineer.  Principal  contributors  to  the  project  and  their  associated 
areas  of  responsibility  Include:  Messrs.  D.  George  and  G.  R.  Schriro  - 

Overall  Program  Integration  and  Final  Checkout;  Dr.  J.  Markowitz  - Integration 
of  Flutter  Redesign  and  Strength  Redesign  Program  Functions;  Messrs.  E.Lerner 
apt!  J.  H.  Berman  - Evaluation  of  Candidate  Flutter  Redesign  Procedures;  Messrs. 

R.  R.  Chipman  and  M.  Chemoff  - Development  of  Integrated  Flutter  Analysis  Module; 
Dr.  W.  J.  Dwyer  - Strength  Analysis  and  Redesign  Module;  Mr.  P.  Shyprykevich  - 
Applied  Loads  Analysis  Module;  Messrs.  M.  J.  Shapiro  and  S.  Goldenberg  - 
Vibration  Analysis  Module.  The  continued  assistance  and  advice  of  Mr  J.  Smedfjeld 
and  Capt.  S.  M.  Batill  have  been  greatly  appreciated.  The  authors  also  wish  to 
acknowledge  Mr.  W.  Mykytow  and  Dr.  L.  Berke  for  initiating  this  effort  and  for 
their  valuable  suggestions  during  the  course  of  the  project. 
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INTRODUCTION 

In  order  to  neet  the  increasingly  restrictive  budgetary  and  schedule 
constraints  imposed  on  aerospace  vehicle  development  programs,  considerable 
attention  has  been  devoted  in  recent  years  to  automation  of  the  aircraft 
design  process.  This  effort  has  largely  focused  on  the  refinement  and  inte- 
gration of  existing  analysis  tools  (e.g.,  see  References  1-1  and  1-2).  At 
the  same  time,  some  significant  advances  have  been  made  in  the  development 
of  automated  redesign  procedures  as  a further  means  of  accelerating  the 
the  design  cycle  (e.g.,  see  Reference  1-3). 

The  computer  program  described  in  this  report  falls  into  both  of  these 
categories.  Its  capabilities  encompass  integrated  interdisciplinary  anal- 
ysis as  well  as  enhanced  automated  redesign  for  aircraft  lifting  - surface 
structures.  The  analysis  capability  includes  loads  prediction,  structural 
analysis,  vibration  mode  determination,  and  flutter  analysis.  The  redesign 
procedure  minimizes  the  weight  of  a lifting  - surface  structure  in  the 
presence  of  combined  strength  and  flutter-speed  constraints.  The  entire 
program  is  known  as  FASTOP,  for  .Flutter  J^nd  STrength  Optimization  Program. 

Before  describing  the  specific  features  of  this  program,  it  is  valuable 
to  review  some  of  the  background  and  objectives  which  governed  its  devel- 
opment . 

The  strength  analysis  and  redesign  module  in  FASTOP  was  developed 
several  years  ago  under  a contract  sponsored  by  the  Air  Force  Flight  Dynamics 
Laboratory.  This  program,  known  as  ASOP  (Automated  Structui al  Optimization 
jhrogr  ~m) , automatically  resizes  the  gages  of  a finite-element  structures 
model  to  achieve  a fully  stressed  (near-minimum-weight)  design  for  spec- 
ified design  load  conditions.  Prior  to  the  advent  of  automated  strength 
resizing  procedures,  the  stress  analyst  was  faced  with  the  time  consuming 
task  of  computing  element  stresses  based  on  the  results  of  a finite-element 
analysis,  and  then  manually  resizing  the  preliminary  gages  of  the  structures 
model.  This  process  bad  to  be  repeated  until  an  acceptable  design  was 
achieved.  Consequently  the  automated  resizing  capability  in  ASOP  has  resulted 
in  a significant  reduction  in  the  time  required  for  the  strength  design  of 
airplane  structures. 
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The  major  objective  in  developing  FASTOP  has  been  to  integrate  the 
aforementioned  strength  redesign  program  with  a new  automated  procedure  for 
minimum-weight  resizing  to  satisfy  a flutter-speed  constraint.  This  require- 
ment originated  from  the  obvious  inefficiencies  in  existing  flutter  preven- 
tion procedures,  wherein  the  flutter  analyst  relies  largely  on  judgement  in 
pursuing  an  adequate  flutter  "fix".  Such  an  approach  often  leads  to  many 
trial  and  error  studies  which  are  particularly  inefficient  because  of  the  non- 
automated  interface  between  the  flutter  analysis  and  structural  analysis 
procedures.  That  is,  any  stiffness  change  to  be  considered  in  the  course  of 
a flutter  investigation  must  be  evaluated  through  manual  changes  to  the  input 
data  of  the  structural  analysis  program.  New  stiffness  properties,  vibration 
modes,  and  a flutter  speed  are  then  computed  on  a step-by-atep  basi3  via 
individual  computer  programs.  The  time  required  for  each  of  these  tasks 
and  the  number  of  "fixes"  to  be  evaluated  obviously  increase  in  proportion 
to  the  complexity  of  the  structures  model  under  investigation. 

Automation  of  this  interactive  strength/flutter  redesign  process  has 
been  accomplished  in  FASTOP,  providing  the  user  with  redesign  capability 
for  small-or  large-scale  structures  models.  Moreover,  the  flutter  redesign 
procedure,  which  evolved  from  an  extensive  evaluation  of  candidate  app- 
roaches, has  been  demonstrated  to  achieve  a minimum-weight  redesign  for 
flutter-critical  configurations. 

This  report  describes  the  various  analysis  and  redesign  procedures 
included  in  FASTOP.  Results  obtained  from  the  analysis  and  resizing 
of  three  sample  structures  are  also  presented.  After  first  providing 
an  overview  of  the  entire  FASTOP  system,  each  of  the  major  analysis 
and  redesign  functions  is  discussed  in  a separate  section.  It  will 
be  noted  that  the  depth  of  detail  varies  in  each  of  these  sections. 

Where  adequate  documentation  of  well-known  methods  already  exists, 
as,  for  example,  in  the  discussion  of  the  structural  analysis  methods, 
considerable  reliance  is  placed  on  this  documentation,  and  the  discussion 
only  summarizes  the  procedures.  On  the  other  hand,  •’ven  though  some 
other  techniques  are  already  documented,  it  is  felt  that  they  may  be 
less  familiar  to  the  reader  and  that  he  would  benefit  from  having  them 
discussed  in  detail  and  included  in  this  volume  for  completeness. 

In  other  instances,  where  the  analysis  capability  provided  in  a parti- 


cular  discipline  extends  into  areas  beyond  the  original  intended  scope  of 
this  contract,  as,  for  example,  in  the  consideration  of  wing-body  aero- 
dynamic interference  and  in  the  computation  of  divergence  speed  in  the 
flutter  analysis  section,  drily  a brief  discussion  of  these  featurc-3  is 
included . 


i 


! 

i 


» 


t 

r 


3 


» 


bs^iaauhLt. 


.n  m rY 


• f-  * > ' t o * - * HK  //'  ‘ v 


Section  2 

OVERVIEW  OF  THE  FASTOP  SYSTEM 

The  FASTOP  system  provides  capability  for  the  analysis  and  near-minimum- 
weight  structural  sizing  of  a lifting  surface  to  meet  strength  and  flutter- 
speed  requirements.  A functional  flow  diagram  of  the  system  is  presented  in 
Figure  2.1.  The  package  is  comprised  of  two  major  programs,  each  one  designed 
to  perform  succe&'.ive  analysis  and  resizing  functions  in  a single  computer 
submission.  The  ji'crength  Optimization  program  (SOP)  focuses  on  basic  aspects 
of  static  structural  analysis  and  minimum-weight  design  for  strength  require- 
ments. It  provides  for  automated  calculation  of  applied  loads,  performance 
of  conventional  strength  and  flexibility  (or  stiffness)  analysis,  and  auto- 
mated resizing  cf  a structural  idealization  to  achieve  a fully  stressed  de- 
sign.. It  also  prepares  data  required  for  direct  input  to  the  second  major 
program.  The  flutter  Optimization  program  (FOP)  addresses  dynamic  analysis 
requirements  and  provides  the  redesign  capability  for  achieving  a desired 
value  of  flutter  speed  with  minimum  cost  in  weight.  Using  output  data  from 
the  first  program,  it  proceeds  to  establish  mass  matrix  input  for  vibration 
mode -analysis,  compute  normal  mode  shapes  and  frequencies,  determine  the 
suri.'oe’s  critical  flutter  speed,  and  perform  resizing  if  desired  to  in- 
crease flutter  speed.  Finally,  the  second  program  saves  data  required  for 
re-entering  SOP. 

The  use  of  this  two-program  approach  for  redesign  proceeds  as  follows. 
First,  using  SOP,  the  structure  is  sized  to  satisfy  its  strength  requirements 
with  a fully  stressed  design.  The  structural  elements  after  this  step  can  be  in 
either  of  two  categories  - fully  stressed  (i.e.,  3trength-critical)  or  at  a 
prespecified  minimum  gage  (as  dictated,  for  example,  by  manufacturing  consid- 
erations). Then,  using  FOP,  resizing  of  structural  and/or  mass-balance 
variables  is  performed  to  increase  the  critical  flutter  speed  of  the  surface. 

No  structural  elements  can  be  reduced  in  this  step,  since  this  would  cause 
them  to  fall  below  the  gage3  previously  required  by  SOP.  The  structural 
variables  that  are  increased  in  the  second  step,  plus  any  mass-balance  variables 
that  have  been  specified  by  the  user,  constitute  an  .initial  set  of  "flutte’.’- 
critical"  elements,  i.e.,  elements  whose  sizes  have  been  dictated  by  flutter- 
speed  requirements.  The  structural  variables  that  have  been  increased  by  FOP 
are  now  removed  from  the  previous  strength-critical  end  minimum-gage  categories. 
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Figure  2.1.  Functional  Flow  Diagram  for  the  FASTOP  System 


Since  the  resizing  of  certain  structural  elements  to  achieve  an  increase 
in  flutter  speed  may  alter  the  internal  load  distributions  and  therefore  the 
resulting  gage  requirements  for  strength  considerations,  SOP  is  now  re- 
entered. This  is  the  first  step  toward  achieving  a minimum-weight  design 
accounting  for  strength/ flutter  interaction.  In  this  step,  no  flutter- critical 
element,  and  of  course  no  element  at  the  prescribed  minimum  gage,  can  be 
resized  downward,  but  the  remaining  elements,  which  constitute  the  newly 
defined  set  of  strength-critical  elements,  can  be  either  decreased  or  increased. 
Again,  a reclassification  of  the  various  elements  into  stfength-critical, 
flutter-critical,  and  minimum-gage  categories  takes  place.  For  example,  if 
a previously  flutter-critical  element  is  resized  upward  to  meet  strength 
requirements,  it  is  now  moved  into  the  strength-critical  set;  similarly,  if 
a previously  3trength-critical  element  can  be  reduced  in  weight  until  it  reaches 
its  prescribed  minimun  gage,  it  is  now  put  into  the  minimum-gage  category. 

At  this  point,  FOP  is  entered  for  the  second  time,  and  the  strength/ 
flutter  interactive  redesign  continues.  The  resizing  that  can  take  place  in 
this  second  pass  through  FOP  is  more  flexible  than  that  which  occurred  the  first 
time,  in  that  now  there  is  a set  of  flutter-critical  structural  elements  that 
can  be  resized  either  upward  or  downward.  Any  downward  resizing  cannot  violate 
the  values  required  by  strength  or  minimum-gage  considerations,  however.  As 
before,  elements  in  the  new  sets  of  strength-critical  and  minimum-gage 
variables  can  only  be  resized  upward,  and  if  this  occurs  they  are  reclassified 
as  flutter-critical. 

Subsequent  interactive  application  of  the  two  programs  proceeds  in  a 
manner  similar  to  the  second  passes  until,  in  the  opinion  of  the  user,  the 
process  is  sufficiently  converged.  The  final  design  will  consist  of  a set  of 
flutter-critical  elements  which  have  nearly  uniform  flutter- velocity  derivatives, 
a set  of  strength-critical  elements  which  are  fully  stressed,  and  a set  of 
elements  which  are  at  the  user  prescribed  minimum  gages. 

While  it  is  felt  that  the  two-program  approach  provides  the  most  logical 
''stop  points"  needed  by  the  user  for  appraisal  of  intermediate  results  and 
possible  revision  of  redesign  step-size  control  parameters,  it  is  still  possible 
to  use  a multiple-step  option  to  exercise  the  analysis  and  redesign  functions 
of  both  programs  without  interruption.  It  is  also  noted  (as  indicated  in 
Figure  2.1)  that  prevision  is  made  to  exit  either  program  after  performing 
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specific  analysis  or  redesign  functions,  thereby  affording  the  user  the  oppor- 
tunity to  monitor  the  results  even  more  closely  or,  if  he  wishes,  exercise 
individual  portions  of  the  system  independently.  Finally,  if  a user  observes 
that  flutter  resizing  occurs  in  a vary  local  area  with  only  minor  interaction 
with  strength  requirements,  the  flutter  resizing  program  can  be  used  several 
times  in  succession,  using  coupled-mode  flutter  analyses,  before  returning  to 
another  strength  resizing  step.  In  this  situation  it  is  recommended,  however, 
that  SOP  be  used  in  an  analysis  mode  to  compute  a new  flexibility  (or  stiffness) 
matrix  for  a new  normal-mode  calculation  after  each  resizing  step;  if  coupled 
modes  are  used,  the  accuracy  of  the  flutter-velocity  derivatives  deteriorates 
rapidly,  and  improper  resizing  steps  can  result. 

In  both  programs,  considerable  emphasis  has  been  placed  on  the  modular 
programming  concept,  so  that  the  system's  capability  can  readily  be  extended 
in  the  future.  The  analysis  and  redesign  functions  of  both  programs  are 
performed  with  the  following  six  modules: 

• Applied  Loads 

• Structural  Analysis  and  Resizing 

• Transformation  Procedures 

• Vibration  Analysis  (including  Mass  Matrix  Definition) 

• Flutter  Analysis 

• Flutter  Resizing 

The  capabilities  of  the  individual  analysis  and  redesign  routines  are  briefly 


summarized  below: 

t 

Applied  Loads 

Aerodynamic 

Maximum  number  of  flight  conditions 

(subsonic  and/or  supersonic)  8 

Maximum  number  of  control  surfaces  8 

‘aximum  number  of  aerodynamic  panels  100 

Inert .al 

Maximum  number  of  flight  conditions  8 

Maximum  number  of  distributed  (point)  masses  1000 

Maximum  number  of  concentrated  (mass  and  inertia)  masses  IOC 


Structural  Analysis  and  Resizing 

Primarily  for  metallic  structures  (limited  composites  capability) 


Maximum  allowable  number  of  finite  elements 

to  define  the  structures  model  3000 

Maximum  number  of  structures  model  node  points  1000 

Maximum  number  of  structures  model  degrees 

of  freedom  6000* 

’•■Maximum  number  of  applied  load  conditions  8 

Vibration  Analysis 

Maximum  number  of  dynamics  model  degrees 

of  freedom  200 

Flutter  Analysis 

Assumed-pressure-function  and  doublet- 


lattice  routines  for  subsonic  flow  M=0-*»>0.9 

Mach-box  routine  for  supersonic  flow  M=1.3-*-3.0 

Maximum  number  of  modes  for  flutter  analysis  20 

Maximum  number  of  control  surfaces  on  main 

surface  for  doublet-lattice  and  Mach-box  routines  5 

Maximum  number  of  aerodynamic  panels: 

Doublet-lattice  4 00 

Mach-box  350 

Flutter  Resizing 

Maximum  number  of  elements  which  can 
be  resized  for  flutter: 

Structural  elements  2000 

Mass-balance  elements  20 


The  following  sections  describe  the  theory  and  procedures 
in  each  program  module.  Each  section  begins  with  a brief  summary  to  enable 
the  reader  to  quickly  grasp  the  intent  of  each  analysis  and  to  understand 
its  relationship  to  the  overall  system. 

♦Reduces  to  3000  if  subsequent  flutter  resizing  uses  a free-free  vibration 
model;  unchanged  for  cantilever  model. 
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Section  3 

.APPLIED  LOADS  ANALYSIS 


3.1  SUMMARY 

The  applied  loads  analysis  module  provides  the  capability  for  computing 
aerodynamic  and  inertial  loads  for  specified  maneuver  conditions.  The  aero- 
dynamic forces  are  computed  by  modeling  the  lifting  surface  with  a distribution 
of  vortices  for  subsonic  flew  or  with  a distribution  of  sources  for  supersonic 
flow.  For  flight  conditions  expressed  in  terms  of  vehicle  load  factors  and 
angular  accelerations  and  velocities,  inertial  loads  are  computed  at  a grid- 
work  of  concentrated  and  distributed  masses.  The  aerodynamic  and  inertial 
forces  are  then  transformed  from  their  respective  math  models  to  the  structures 
model  using  transformations  computed  in  a separate  transformation  analysis 
module  described  in  Section  5.  FASTOP  also  provides  the  capability  for  the 
direct  input  of  known  loads  in  the  structures  model. 

3.2  AERODYNAMIC  LOADS 

To  obtain  surface  aerodynamic  forces,  the  planfonn  is  subdivided  into 
an  arbitrary  number  of  small  trapezoidal  panels  (not  more  than  100 ) in  a 
fashion  dictated  by  the  overall  planform  geometry  and  the  locations  of  the 
control  surfaces.  The  number  of  panels  in  the  chcrdwise  direction  can  vary 
over  the  span.  Using  the  same  panel  geometry  for  all  Mach  numbers,  aero- 
dynamic influence  coefficients  corresponding  to  these  panels  are  computed  and 
stored  on  tape  using  either  subsonic  vortex-lattice  theory  (Reference  3”l)  or 
the  supersonic  source  distribution  theory  (Reference  3-2).  For  subsonic  flow 
only,  the  effect  of  a fuselage  on  the  lift  on  the  wing  can  be  modelled  in 
this  program  using  the  method  of  images.  Ry  multiplying  the  aerodynamic 
influence  coefficients  by  prescribed  dynamic  pressures  and  downwash  distribu- 
tions, the  forces  are  determined  for  the  selected  loading  conditions.  Pro- 
vision for  various  correction  factors  is  included. 


3.2.1  Subsonic  Influence  Coefficients 

In  subsonic  flow,  the  calculation  of  aerodynamic  influence  coefficients 
is  based  on  the  vostex-lattice  method  of  Reference  3-1.  If  a distribution  of 
vortices  is  placed  on  the  surface  of  a planform,  the  resulting  downwash  at 
iny  point,  P^,  i3  related  to  the  circulations  of  the  vortices  by: 

w (V " § " life  p)  K(?j»  p» M>  ds>  o.u 

*S 

where 

P is  any  point  on  the  planform  with  coordinates  5,  I),  C 
Pj  is  a point  with  coordinates  x,  y,  z 

x,  § are  streamwise  coordinates 

y,  T1  are  spanwise  coordinates 

z,  Q are  vertical  coordinates 

S is  the  surface 

w is  the  downwash  angle 

w is  the  downwash  velocity 

U is  the  free  stream  velocity 

Y is  the  circulation 

M is  the  Mach  number 

K is  the  kernel  function  representing  the  downwash  created  at  a point 
due  to  a unit  circulation  over  a unit  area  of  a vortex  sheet. 

From  the  Kutta-Joukowsky  theorem,  it  is  known  that  circulation  can  be  related 
to  the  lift  and,  hence,  to  the  pressure  coefficient,  C^.  Consequently,  the 
above  equation  can  be  written 

w (V  = A JfC p(P)  K(Pj’  p>  M)  dS-  (3>2) 
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In  the  present  method,  the  given  surface  is  divided  into  a lattice  of  I 
panels,  the  sides  of  which  parallel  the  free  stream  (see  Figure  3.1)  and  over 
each  of  which  the  pressure  is  assumed  constant.  Equation  (3.2)  becomes 


K(P  , P,  M)  dS. 


Panel  i 


With  the  further  assumption  that  the  vortices  in  each  panel  are  condensed  to 
a single  horseshoe  vortex,  the  bound  portion  of  which  li js  at  the  quarter  chord 
of  the  panel,  the  equation  reduces  to 

* (Fj)  - gfc  £ % ASi  ftlVy  ? (5l),  M)  dS,  (3.4) 

i=1  1 "Panel  i 

where  a§  is  the  average  chord  length  of  the  i*'*1  panel  and  the  integration  is 
taken  along  the  quarter-chord  of  the  i panel.  Using  the  Biot-Savart  law, 
the  final  integral  is  evaluat-  d in  Reference  3-3,  yielding 

/K<V  P (5i>-  *>  as  " e31  * FJi  ♦ V (3.5) 


Panel  i 


where 


Eji  is  the  downwash  at  point  j induced  by  a horseshoe  vortex  of  unit 
strength  and  length  at  panel  i 

F^  is  the  contribution  to  due  to  the  trailing  vortices 

is  the  contribution  to  due  to  the  bound  vortex. 

If  it  is  desired  in  determining  the  pressure  distribution  on  a wing  to 
account  for  the  effects  of  the  presence  of  a fuselage,  this  is  accomplished  as 
in  Reference  3-4  by  including  within  the  fuselage  an  image  of  each  horseshoe 
vortex  on  the  wing  (see  Figure  3.2).  In  this  idealization,  the  fuselage  is 
assumed  to  be  a circular  cylinder  of  radius  R and  the  image  of  a point  (x,  y) 
is  located  at  (x,  R ,/y).  The  downwash  induced  by  the  point  and  its  image  is 


2 = ('F 

Si  K Ji 


* F5i>  * <c3i  * °ji> 


(3.6) 


Ji 


V Z*  '/!  •//  V/  •//  V7  7 I'll  'I • 


where  F'  and  G*  are  due  to  the  linage  and  where  the  correction  factor, 
J*  0 * 

r\  g 

(l  + K7y*),  on  the  dowr.wash  due  to  the  bound  vortex  is  derived  from  two- 
dimensioned  theory. 


I 


Figure  3*2.  Cylinder  with  a Fair  of  Horseshoe  Vortices 
and  their  linages. 
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Substituting  the  appropriate  value  of  E^i 
one  obtains 


for  the  integral  in  Equation 


i 

8n 


I 

E 

i=l 


AS,  e, 


i di* 


(3.7) 


If  a point  is  chosen  at  the  3A-chord  and  midspan  of  each  of  the  I panels  and 
the  above  equation  is  applied  at  each  such  point,  a system  of  simultaneous 
linear  equations  results,  which  in  matrix  notation  is 


M - A M H {si- 


This  system  can  be  solved  for  the  pressure  distribution: 

l-l 


H- 


(3-9) 


From  this  equation,  an  aerodynamic  influence  coefficient  matrix  can  be 
defined  as 


(3.10) 


The  downwash  angle  distribution,  fw},  used  in  Equation  (3.9)  can  be  con- 
sidered to  consist  of  contributions,  (wj,  due  to  the  rigid  surfaces  inclina- 
tion to  the  free  stream  and  contributions,  Iwg},  arising  from  the  deflections 
of  any  control  surfaces  on  the  surface.  The  first  of  these  contributions  is 
comprised  of  three  terms 

{vjJ  » {«,}  + *2  («2i  + k3  (®3)>  (3.11) 


where 

(or  } is  the  effective  angle  of  attack  of  the  surface  optionally  in  the 
presence  of  a fuselage 

{»2}  is  the  distribution  of  local  incremental  angles  of  attack  due  to 
camber  and  twist 


{cr^}  is  a distribution  of  additive  corrections  to  the  local  incidences 
based  on  empirical  or  experimental  data 

kg,  k^  are  scalar  correction  factors  also  based  on  empirical  or  experi- 
mental *__ta. 

The  effective  angle  of  attack  equals  the  sum  of  tne  geometric  angle  of  attack 
of  the  wing  relative  to  the  fuselage  and  terms  due  to  the  fuselage's  inclina- 
tion and  the  upwash  induced  by  this  inclination. 

a 

t*e)  to  + af  [1  + FfVy2]  [cx  J {1},  (3.12) 

where 


is  the  geometric  angle  of  incidence  of  the  wing  relative  to  the 
fuselage  reference  line 

ofj,  i3  the  angle  of  attack  of  the  fuselage  reference  lire 
R is  the  mean  radius  of  the  fuselage 

y is  the  spanwise  coordinate  of  the  panel  of  interest  measured  from 
the  fuselage  centerline 


C,  1 are  scalar  corrections  factors. 


The  contribution  to  the  downwash  due  to  N control  surface  deflections  is 

^ ' n?i(Sl*n  tC®3n  [U)”  * \ C#2>ni) 

*I(v  n„  (v}). 


(3.13) 


where 


N 


is  the  number  of  control  surfaces 

th 


5 ,6  are  the  rotations  of  the  n left  and  right  control  surfaces 

L>n  «»n  respectively 

til 

(U}  are  fractions  (OsUj,n  s.  l)  denoting  the  portion  of  the  j 


til 

aerodynamic  panel  that  lies  on  the  n control  surface, 

f } are  additive  corrections  to  the  locsl  rotation  of  the  n 
* <}n 

control  surface 


th 


k^  are  scalar  correction  factors, 
n 
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It  should  be  noted  that,  in  the  absence  of  any  correction  factor,  the 
dovmvash  distribution  on  a wing  without  fuselage  or  control  surfaces  is  quite 
simple,  namely  {w}  = or^  {l}  + {«2}. 

3.2.2  Supersonic  Influence  Coefficients 

In  supersonic  flow,  aerodynamic  influence  coefficients  are  obtained  as 
described  in  Reference  3-2  using ‘a  distribution  of  sources.  The  pressure  at  a 
point  on  a planar  surface  is  related  to  the  velocity  potential  and,  thence, 
to  the  downwash  by 

W • 1 ■ ? iffi  4»"-  (3-w 

s 


Q 

p is  the  pressure  coefficient 

Pj  is  a point  with  coordinates  x and  y 

x and  5 ore  streamwise  coordinates 
y and  T|  are  spar.wise  coordinates 
U is  the  free-stream  velocity 

cp  is  the  velocity  potential 

w is  the  downwash  angle  at  5,  1} 

R « V(x  - 5)2  - p2  (y  - t)2 
P2  « & - 1 

M is  the  Mach  number 

S is  the  area  of  integration  bounded  by  the  inverse  Mach  cone 

emanating  from  (x,  y). 

For  a surface  with  supersonic  leading  and  trailing  edges,  the  area  of 
integration  is  bounded  by  the  leading  edge  (see  Figure  3- 3a).  Near  a side 
edge  of  such  a surface,  however,  the  potential  is  influenced  by  a region  off 
the  wing  as  indicated  by  SD  in  Figure  3.3b.  For  such  a case,  Reference  3-5 
shows  that  the  effect  of  this  off -wing  region  is  to  negate  +he  conti tbution  of 
the  on-wing  region  delimited  by  the  leading  edge,  the  side  edge  and  the 
reflected  Mach  line.  Consequently,  only  the  shaded  area  Sw  need  be  considered 
for  the  integration  in  Equation  3.13. 


| This  concept  can  be  extended  to  the  case  of  subsonic  leading  edges, 

j Referring  to  Figure  3* *+»  the  areas  between  the  foremost  Mach  waves,  OB'  and 

i OD'  and  the  wing  leading  edges,  OEB  and  0(2),  will  affect  the  potential  and 

! cancel  the  influence  of  the  regions  EEC,  EFH,  etc.  ahead  of  the  reflected 

Mach  lines  BC,  EF,  etc.  Thus,  the  areas  to  be  considered  in  the  integration 
are  the  successively  forward  quadrilaterals  A BCD,  CEFG,  etc. 

i 

t 

i 


Figure  3.*+.  Integration  Areas  for  a Ming  with 
Sub, sonic  Leading  Edges. 


Returning  to-  Equation  (3,Xlf } and  performing  the  differentiation  indicated, 
Ot»j obtains" 

(3-15) 

J yS'  i-.E. 

which  can  he  simplified  to 

W-liTt  TS*s/id,>-  <3-l6> 

S L.E. 


In  the  above  equations,  T.  E.  and  L.  E.  signify  that  the  line  integral  is  to 
be  evaluated  along  the  trailing  and  leading  edges  of  the  area,  S,  respectively. 

To  evaluate  the  above  equation,  the  planfona  is  divided  into  a number  of 
trapezoidal  panels,  as  was  done  in  the  subsonic  case.  Over  each  panel,  the 
downwash  is  assumed  to  be  constant.  Consider  a typical  panel  and  its  neighbors 
as  depicted  in  Figure  3.5.  The  contribution  to  the  pressure  at  (x,  y)  of  the 
wedge  aft  of  A2C  is 


ACG 


2 CC  6w  dgdT)  2 (f  w 

v JJ  as  R n / r 


dll 


ACG 


ABC 


mS*f  -11*11  *f. 

ACG  ABC  ACFD  DFG  ABC 


(3-17) 


The  contribution  from  the  wedge  aft  of  DEF  is 


(3.18) 


Subtracting  these  two  equations,  one  obtains  the  contribution  of  the  strip  ACFD: 

hem-ff 


(3.19) 


Figure  3*5.  Planform  Divided  into  Panels. 

This  can  be  broken  up  into  the  contributions  associated  with  panel  i (shaded  in 
Figure  3.5)  and  those  associated  with  its  spanwise  neighbor(s).  For  panel  i, 
one  has: 


•//*/■/■#//  ’ <3'20) 
ABED  AB  DE  Si  L.E^  T.E^ 

where  only  those  portions  of  S^,  L.E.^  and  T.E.^  that  are  within  the  Mstch  cone 
are  considered.  Since  the  downwash  is  assumed  constant  over  each  panel,  the 
surface  integral  in  Equation  (3.20)  equals  zero  and  the  net  result  involves 
only  line  integrals  along  the  leading  and  trailing  edges.  Hence,  Equation 
(3.16)  becomes: 
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(3.21) 


where  a Is  the  number  of  panels  within  the  inverse  Mach  cone.  The  regaining 
integral sare  easily  evaluated,  since  the  equations  for  the  panel  leading  and 
trailing  edges  are  those  of  straight  lines. 

For  thin  airfoils,  the  magnitude  of  the  dovnwash  on  both  sides  of  the 
surface  can  be  assumed  equal  to  the  local  angle  of  attack  of  the  mean  surface. 
Hence,  the  net  pressure  coefficient  is,  merely,  twice  the  value  given  in  Equa- 
tion (3.21)  with  the  dovnwash  evaluated  using  Equation  (3.11).  Fuselage  effectc 
are  not  included,  however,  since  the  formulation  given  previously  is  not 
adequate  for  supersonic  flow. 

Denoting  the  two  line  Integrals  in  Equation  (3.2l)  as  A^,  the  net  pressure 
becomes 


VV  -sin 


i=m+l 


(3.22) 


where  n is  the  total  number  of  surface  panels  and  the  second  summation  (equal 
to  zero)  contains  those  panels  not  within  the  inverse  Mach  cone  emanating  from 
P^.  Applying  Equation  (3*22)  to  the  center  of  each  panel,  one  obtains 

{cp}*  [aic]  H'  (3,23) 

where 

[aiC]  « | [a]  . (3.24) 

The  downvash  distribution  to  be  used  in  Equation  (3.23)  is  specified  in  accor- 
dance with  Equations  (3.H),  (3.12),  and  (3.13).  It  is  noted  that  the  fuse- 
lage radius,  R,  must  be  entered  as  zero  for  supersonic  loads  analysis. 
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3.2.3  Aerodynamic  Forces 


For  & g-*  metrically  synsnetric  planform,  pressure  coefficients  are  computed 
for  the  left  alf  only,  based  on  an  appropriate  combination  of  symmetric  and 
antisynmetric  aerodynamic  influence  coefficients  and  the  specification  of  the 
total  ( genera :iy  asymmetrical)  downvash  angle  distribution  in  terms  of  sym- 
metrical and  antisymmetrical  components: 

ic  i « [aic] 
i p /L  L 

where 

Kym.  + 

Manti-**  [Ml’Mr] 
sym. 

and  L and  R denote  the  left  and  the  right  halves,  respectively,  of  the 
symmetric  plar.form.  With  the  pressure  coefficients  determined,  the  total 
aerodynamic  force  on  each  panel  is  simply 

|cr|  (3.26) 

where  p is  tne  air  density  and  £sj  is  a diagonal  matrix  of  panel  areas. 

Using  the  procedures  described  in  Section  5,  the  aerodynamic  forces  are 
finally  transformed  from  the  aerodynamics  model  to  the  structures  model. 


j w}  + 
* ’sym. 


sym. 


anti- 

sym. 


(3.25) 
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3.3  INERTIAL  LOADS 

Using  D'Alembert's  principle,  rigid-body  inertial  forces  and  moments  at  a 
point  on  a body  can  be  written  as 

..  • 

-F  « m r + mwx  'jjxpJ+miuxp 

and  -M  * H + u)  x H,  (3.2' 

where 

F is  the  inertial  force  vector 
M is  the  inertial,  moment  vector 

r is  the  position  vector  locating  the  origin  of  the  reference  axes 
on  the  body 

u is  the  mass  of  the  point 

(i>  is  the  angular  velocity  of  the  reference  axes 
p is  the  position  of  the  point  relative  to  the  reference  axes 
H is  the  moment  of  momentum  of  the  point. 

Expansion  of  these  equations  yields 

FX,  i ■ (xt  - XQ)  (Q,2  + R2)  m^  + (yt  - yQ)  (R  - PQ/  ^ 

-(^  - zQ)  (RP  + 4 ) • g Sj 

Fy,  i * “ y0}  (r2  + + (zi  - z0)  mi 

e 

-(x.  - x ) (PQ  + R)  m.  - g N m. 


I.  1 


n . 

y*  * 


n . 

z,  i 


* (2i  ■ zo)  ^ + Q2)  mA  + (xA  - x0)  (i  - RP)  mA 
“(yA  - yd)  (OR  + P)  aA  - g Nz  nA 

• -Jxx,  i * - Xxy,  i <»  “ * Xxz,  i ^ + « 

1 ^ ' “*>  - (I«,  i - Jyy,  i>  «»  . 

■ _Iyy.  i « - t (or  - « * , «*  * p) 

•1b.  i ^ ' "*>  ' »«,  i - i>  w 

■ -1..,  i * - W <*>  - » * i » * « 
-V,  i ' t?>  - (1n-,  i ' 'x,,  i>  «>  • 


(3.28) 


where 

g is  the  acceleration  of  gravity 

*i»  y1#  Zji  aA  are  the  coordinates  ond  mass  of  the  ith  point 

Ixx,  i*  Iyy<  i*  etc.  are  the  inertia  properties  about  the  center  of 
gravity  of  the  ith  point 

N , W » N are  load  factors 
a 4 

P,  Q,  R are  x,  y,  z components  of  the  angular  velocity,  ui,  respectively 

XQ>  y0>  z0  are  the  coordinates  of  the  reference  point  (usually  the  center 
of  gravity  of  the  airplane). 


There  are  two  types  of  mass  data  that  the  user  may  provide  for  inertial 
load  computations.  The  first  type,  termed  "distributed  masses",  are  simply 
point  masses  with  no  rotational  inertia  properties.  Thus  the  output  of  the 
inertial  loads  routine  for  distributed  masses  consists  only  of  forces  st  the 
center  of  gravity  of  each  mass  item.  Distributed  masses  are  most  typically 
used  to  define  the  mass  of  the  finite  element  structures  model  (including 
non-optimum  factors  where  appropriate)  and  are  logically  assigned  to  struc- 
tures model  node  points.  The  inertia  properties  of  the  structure  are  im- 
plicit in  the  node  point  mass  latributlon,  The  second  type  of  mass  data, 
referred  to  as  "concentrated  masses",  have  both  mass  and  inertia  properties. 
In  this  case  the  inertial  loads  routine  computes  inertial  forces  and  moments 
at  the  center  of  gravity  of  each  mass  item.  Concentrated  masses  are  used 
to  represent  large  masses  such  as  engines,  stores,  etc. 

The  computed  inertial  loads  acting  at  the  centers  of  gravity  of  the 
various  masses  are  subsequently  distributed  to  the  n.jde  points  of  th';, 
structures  model  using  the  transformation  procedures  dencribed  in  Section  5. 

3.4  PROVISION  FOR  DIRECT  INPUT  OF  APPLIED  LOADS 

For  non-maneuver  loading  conditions  (e.g. , landing  or  gust  loads)  or 
maneuver  loading  conditions  for  which  applied  loads  have  been  determined 
directly  from  test  data,  separate  loading  conditions  consisting  of  forces 
at  structural  nodes  may  be  prescribed  by  the  user. 
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Section  4. 

STRUCTURAL  ANALYSIS 


4.1  SUMMARY 

Structural  analysis  and  strength  resizing  capability  are  provided  by 
utilizing  the  Automated  Structural  Optimization  Program  (ASOP)  described  in 
detail  in  Reference  4-1.  The  analysis  procedure  employs  the  matrix  displace . 
ment  method  to  obtain  nodal  deflections  and  internal  element  comer  forces . 

The  comer  forces  are  then  converted  into  forces  that  are  equivalent  to  those 
obtained  by  the  matrix  force  method,  and  these  are  subsequently  transformed 
into  representative  element  stress  levels.  Resizing  for  strength  requirements, 
in  the  presence  of  minimum-gage  constraints,  is  accomplished  through  the  itera- 
tive use  of  ratios  of  actual-to-allowable  stress  to  satisfy  the  fully-stressed- 
design  optimality  criterion. 

The  following  paragraphs  provide  a summary  of  the  finite  elements  that 
can  be  used  in  a structural  idealization  within  FASTOP.  The  basic  aspects  of 
the  overall  structural  analysis  procedure  are  briefly  reviewed  and  the  numeri- 
cal method  for  solving  the  structure's  load/displacement  equations  is  presented. 
Further  discussion  of  the  redesign  aspects  of  ASOP,  in  the  context  of  the  com- 
bined atrength-and-flutter  optimization  problem,  is  given  in  Section  10. 

4.2  SUMMARY  OF  AVAILABLE  FINITE  EIEMENTS 

The  finite  elements  available  for  structural  modeling  are  those  that  are 
commonly  used  in  major  structural  analysis  involving  metallic  construction. 

With  the  exception  of  the  plate  elements  (i.e.,  the  bending  triangle  and  the 
bending  quadrilateral),  all  of  the  elements  in  the  following  list  have  stiffness 
characteristics  that  are  directly  proportional  to  their  weight,  thus  making 
them  available  for  resizing  in  FASTOP.  Plate  elements  should  be  used  for  anal- 
ysis only.  The  elements,  which  are  discussed  in  Reference  4-1,  include 

• a uniform  bar  element 

a a uniform  beam  element  having  constant  radii  of  gyration  in  its  cross- 
section 

• an  anisotropic  triangular  membrane  element 
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• an  anisotropic  planar  quadrilateral  membrane  element, 
a a quadrilateral  warped  or  planar  shear  panel  element 

• an  anisotropic  warped  quadrilateral  membrane  element 

• a hinged  beam  element 

a an  anisotropic  bending  triangle  element 
a an  anisotropic  bending  quadrilateral  element 
a a combined  triangular  membrane  and  bending  triangle  element 
a a combined  planar  quadrilateral  membrane  and  bending  quadrilateral. 

4.3  REVIEW  OF  ANALYSIS  IROCEDURE 

In  applying  the  matrix  displacement  method,  the  analyst  first  establishes 
a structural  idealization  comprised  of  the  above  elements,  and  representing, 
as  closely  as  possible,  the  actual  topological  arrangement  of  the  primary 
structural  members.  The  required  input  data  is  then  classified  into  groups 
defined  as:  nodal  geometry,  member  data,  boundary  conditions,  material  tables 
and  applied  loads.  The  member  data  contains  both  topological  data  plus  member 
sizes.  Using  this  information,  the  program  assembles  the  total  stiffness 
matrix  [K]  by  superimposing  the  element  stiffness  matrices  compatible  with  a 
global  coordinate  system.  That  is, 


[Kj 


■L-i 

i=l 


(4.1) 


where  a.  is  a function  of  the  design  parameter  t,  and  [k  ] is  the  expanded 
element  stiffness  matrix  for  a unit  value  of  that  design  parameter  with  the 
appropriate  boundary  conditions  applied  to  it.  The  matrices  are  calculated 
only  once  and  stored.  In  succeeding  redesign  cycles  they  are  multiplied  by  the 
n'eir  design  parameter  and  reassembled  to  form  the  new  [K]  matrix. 

Applied  forces  for  the  various  loading  conditions  are  entered  by  referring 
them  to  the  node  point  at  which  they  are  applied  and  to  the  global  coordinate 
direction  in  which  they  act.  The  forces  are  then  transformed  to  correspond  to 
a degree-of- freedom  numbering  scheme  associated  with  the  free  degrees  of  free- 
dom in  the  idealization.  The  resulting  load  matrix,  TP],  is  then  one  in  which 
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the  rows  correspond  to  the  degrees  of  freedom  and  the  columns  correspond  to 
loading  conditions.  The  equations  of  nodal  equilibrium  are  thus 

[K3[d]  « [P],  ' (4.2) 

where  [A]  is  the  matrix  of  nodal  displacements.  This  system  is  solved  by 
using  a modified  version  of  the  Cholesky  algorithm  (discussed  in  the  next  sub- 
section) and  the  resulting  displacements  are  then  converted  into  element  comer 
forces,  [q],  by 


fq3  * CS^tA],  (4.3) 

in  which  [S^]  embodies  both  element  stiffness  matrices  and  appropriate  force 
transformations  from  a global  to  a local  element  coordinate  system. 

The  procedure  adopted  for  defining  internal  stress  levels  is  the  "nodal 
stress  method"  described  in  References  4-1  and  4-2.  This  procedure  first  con- 
verts the  element  corner  forces,  [q],  into  a new  system  of  forces  that  are 
equivalent  to  those  obtained  by  the  matrix  force  method  of  analysis  (i.e.,  cap 
loads  and  shear  flows).  Then,  an  approximate  strain-compatibility  relationship 
is  used  at  each  structural  node  point  to  determine  the  states  of  stress  at  the 
comers  of  each  finite  element.  Representative  single  values  of  stress  are 
then  computed  at  each  comer  for  use  in  a stress-ratio  resizing  formula  that 
relates  this  value  of  stress  to  its  allowable. 

4.3.1  Solution  of  the  Ioad/Displacement  Equations 

The  solution  of  Equation  (4.2)  is  accomplished  by  employing  the  Cholesky 
algorithm  for  decomposition  of  positive-definite  symmetric  matrices  (see 
Reference  4-3).  This  technique  is  also  used  elsewhere  in  the  FASTOP  system 
(see  Sections  7,8).  The  overall  solution  procedure  involves  the  following 
steps: 

First,  factor  [K]  by  the  Cholesky  algorithm  such  that 

[K]  •=  [L][L]T,  (4.4) 

where  [L]  is  a lower  ^riangilar  matrix.  The  procedure  for  obtaining  [L]  is 
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based  ujr  . 3v  defining  relationship.  Equation  (4.4),  and  the  positive-definite 
symmetric  nature  of  the  stiffness  \trix  [K],  It  first  involves  the  calculation 
of  elements  in  the  first  columr,  ,lven  by 
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where  the  J's  and  k's  are  the  elements  o?'  Cl]  and  [K],  respectively.  There- 
after, each  successive  column  of  [L]  is  computed  with  the  following  recursion 
formulae: 
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Equation  (4.2)  may  now  be  rewritten  as 


[l]Cl3aCa3  - [p3. 


In  the  second  step,  define  [Z3  * C^]  id]  and  solve 


(4.5) 


as 
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[L][Z3  - Cp3 


(4.6) 


for  [Z]  by  successive  substitution,  starting  with  the  first  equation  end  pro- 
ceeding downward. 


For  the  final  step,  solve 


[l3t[a3  * [z] 


for  the  displacements,  [d],  by  successive  substitution,  starting  with  the  last 
equation  and  proceeding  upward. 
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Section  5 


TRANSFORMATIONS  BETWEEN  MATHEMATICAL  MODELS 

5.1  SUMMARY 

Became  different  mathematical  models  are  employed  in  performing  the  major 
analysis  tasks  required  for  the  design  of  a lifting  surface,  it  is  necessary,  in 
an  Integrated  system,  to  provide  automated  methods  that  transfoim  data  obtained 
with  one  model  into  a suitable  form  for  use  in  another.  The  transformation  of 
aerodynamic  and  inertial  applied  loads  into  structural  node  forces  is  achieved 
by  using  fully  automated  "load-beaming"  procedures  that  rely  on  assumed  load 
paths.  These  same  procedures,  with  some  added  special  capabilities,  also  pro- 
vide a basis  for  establishing  a method  that  transforms  flexibility  and  mass  data, 
defined  in  terms  of  e structures  model,  into  a form  for  use  in  vibration  analysis. 

5.2  DEFINITION  OF  REQUIRED  TRANSFORMATIONS 

To  perform  the  analyses  required  for  the  structural  design  of  a lifting 
surface,  four  distinctly  different  mathematical  models  are  usually  employed. 

For  determining  applied  aerodynamic  loads,  an  aerodynamics  model  provides  a 
regular  planar  distribution  of  points  at  which  angles  of  attack  are  specified 
and  forces  normal  to  the  plane  are  calculated.  When  rigid-body  inertial  loads 
are  being  considered,  a weights  model  is  defined,  comprised  of  lumped-mass  points 
representing  both  structural  and  non structural  (i.e.,  leading-  and  trailing-edge 
assemblies,  fuel,  actuators,  etc.)  mass  items.  A structures  model,  used  to  de- 
fine internal  load  distributions  and  elastic  flexibility,  consists  of  an  assem- 
blage of  finite  elements  representing  the  actual  arrangement  of  primary  struc- 
tural members.  To  determine  the  surface's  modes  of  vibration  for  subsequent 
flutter  analysis,  a dynamics  model  is  defined.  This  model  consists  of  lumped 
masses  having  degrees  of  freedom  that  ere  important  for  the  accurate  determi- 
nation of  the  lowest-frequency  vibration-mode  characteristics  and  will  usually 
differ  from  the  weights  and  structures  models. 

It  is  evident  that  an  integrated  analysis  system  must  provide  the  capability 
for  automatically  transferring  data  from  one  model  to  the  next  in  the  analysis 
sequence.  Since  each  analysis  function,  however,  utilizes  a model  tailored  to 
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its  own  specific  discipline,  it  is  also  necessary  to  provide  certain  data  trans- 
formations. The  automated  procedures  included  in  FASTOP  pro'.ide  for 

• transformation  of  loads  in  the  aerodynamics  and  weights  models  to 
equivalent  applied  loadn  in  the  structures  model,  and 

a transformation  of  flexibility  and  mass  matrices  associated  with  the 
structures  model  into  similar  matrices  that  are  representative  of 
the  dynamics  model. 

All  of  these  transformation  are  based  on  principles  of  "load  beaming", 
that  is,  assumed  relationships  that  transfer  loads  in  one  model  to  sets  of 
statically  equivalent  loads  in  another. 

5.3  DESCRIPTION  OF  BEAMING  PROCEDURES 

The  methods  for  transferring  (or  "beaming")  forces  in  the  aerodynamics 
model  and  forces  and  moments  in  the  weights  model  to  the  structural  node  points 
are  illustrated  in  Figures  5.1  and  5.2.  TVo  basic  types  of  beaming  are  provided. 
The  "eight-point"  procedure  is  designed  primarily  to  transfer  applied  aerody- 
namic or  inertial  loads  that  act  at  locations  within  the  geometric  boundaries 
of  the  primary  structure.  The  "four-point"  procedure  is  intended  to  transfer 
loads  that  act  outside  the  structure,  such  as  at  a wing's  trailing  or  leading 
edge.  When  using  either  procedure,  transformation  matrices  are  established 
that  express  loads  at  the  stnntures-model  node  points  (in  its  coordinate  system) 
in  terms  of  unit  applied  loads  in  either  the  aerodynamics  or  the  weights  models 
(in  their  respective  coordinate  systems).  fOr  aerodynamic  forces,  provision  is 
made  for  transferring  only  forces  that  act  normal  to  a reference  plane;  however, 
for  inertial  loads,  all  six  components  of  force  and  moment  may  be  transferred. 

The  program  requires  as  input  the  nodal  geometries  of  the  pertinent  models  and 
correspondence  tables  that  indicate  the  manner  of  beaming  and  the  .nodes  from  and 
to  which  the  loads  are  to  be  transferred. 

The  assumptions  made  in  both  procedures  are  described  in  the  next  two  sub- 
sections. These  are  followed  by  a discussion  of  the  use  of  load  beaming  as  a 
first  step  in  developing  the  transformations  that  are  required  to  convert  data 
from  the  structures  model  to  the  dynamics  model.  It  should  be  noted  in  the 
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following  discussions  that  the  x-y  coordinate  planet  of  all  mathematical  model* 
are  assumed  parallel,  even  though  the  coordinate  systems  may  have  different 
origins, 

5.3.1  Eight- feint  Beaming 

Referring  to  Figure  5.1(a)  or  (b),  point  M represents  an  applied  aero- 
dynamic or  Inertial  load  point;  points  I,  J,  K,  L and  N,  0,  P,  Q represent 
the  upper  and  lower  cover  structural  nodes,  respectively,  to  which  the  applied 
loads  are  being  transferred.  The  "beaming  plane"  contains  point  M and  is  paral- 
lel to  the  structures  model  x-y  plane.  The  line  e-f  is  defined  perpendicular 
to  the  line  a-b,  and  g-h  is  parallel  to  a-b.  Thus,  e-f,  g-h,  and  a normal  to  the 
beaming  plane  fora  a local,  orthogonal  coordinate  system. 

Unit  applied  forces  and  moments  are  first  transformed  into  components  in 

this  local  system.  local  loads  P , P , P , and  M are  first  transferred  along 

« w * y 

line  e-f  to  points  e and  f,  under  the  assumption  that  e-f  acts  as  simple  (pin- 
ended)  beam.  Then,  using  a-b  and  c-d  as  simple  beams,  the  loads  are  transferred 
to  points  a,  b,  c,  and  d,  where  they  are  finally  beamed  to  the  structural  node 
points.  The  local  moments  and  follow  a different  path;  they  are  first 
beamed  to  g and  h,  and  thence  to  a,  b,  c,  and  d by  using  a-c  and  b-d  as  simple 
beams. 

It  should  be  noted  that  all  applied  moments  are  initially  transformed  into 
force  couples  in  the  first  beaming  step,  and  only  force  components  are  eventually 
transferred  to  the  structural  node  points.  Also,  forces  that  are  applied  parallel 
to  a beam  member  are  distributed  to  the  end  points  in  inverse  proportion  to  the 
distances  from  the  point  of  application. 

After  the  forces  at  the  structural  node  points  are  determined,  they  are 
then  rotated  into  the  structural  coordinate  system. 

5.3.2  Pour- feint  Beaming 

As  indicated  previously,  this  procedure  is  particularly  applicable  to  aero- 
dynamic or  inertial  loads  that  are  applied  at  points  external  to  the  structural 
idealization.  Referring  to  Figure  5.2(a),  point  M represents  the  applied  aero- 
dynamic or  inertial  loading  point;  points  I,  J,  K and  L are  the  structural  node 
points  to  which  the  applied  loads  are  being  transferred.  The  "reference  plane" 


is  parallel  to  lines  connecting  I with  L and  K with  J,  and  contains  points  e 
and  f which  are  the  mid-points  of  I-K  and  J-L,  respectively.  Points  a,  b,  c, 
and  d are  orthogonal  projections  of  the  structural  node  points  on  the  refer- 
ence plane.  The  line  g-M  is  perpendicular  to  line  e-f;  these  two  lines  and  a 
common  normal  define  a local,  orthogonal  coordinate  system. 

Unit  applied  forces  and  moments  at  M are  first  transformed  into  this  local 
coordinate  system,  and  then  transferred  to  point  g,  under  the  assumption  that 
M-g  acta  as  a beam  cantilevered  from  point  g.  Next,  they  are  transferred  along 
member  e-f  to  its  end  points,  with  this  member  acting  as  a beam  capable  of 
resisting  bending,  axial  load,  and  torsion.  The  torsion-resisting  capability 
at  e and  f is  assumed  to  be  confined  to  the  planes  connecting  I-K-a-c  and  J-L- 
b-d,  respectively.  The  three  force  components  and  the  concentrated  moment  at 
each  point,  e and  f,  are  then  transferred  to  the  structural  nodes  by  assuming 
that  members  I-K  and  J-L  are  pin-connected  at  their  respective  structural  node 
points.  Where  forces  and  moments  are  applied  parallel  to  a member,  they  are 
distributed  to  the  end  points  in  inverse  proportion  to  the  distances  from  the 
point  of  application. 

For  the  special  ease  where  point  g lies  outside  of  points  e and  f,  as  shown 
in  Figure  5.2(b),  the  torsional  moment  about  e-f  and  the  axial  force  acting  along 
this  member  are  assumed  to  be  resisted  totally  by  the  more  adjacent  support 
point  which  is  f in  the  illustrative  example. 

As  with  the  previous  beaming  procedure,  after  the  forces  at  the  structural 
node  points  are  determined,  they  are  then  rotated  into  the  structural  coordinate 
system. 

5.3.3  Special  Beaming  Features  for  Use  in  Defining  a Dynamics  Model 

lne  preceding  beaming  procedures  are  also  used  in  the  transformation  of 
flexibility  and  mass  data  associated  vith  the  structures  model  into  a form  com- 
patible with  a dynamics  model.  Specifically,  these  procedures,  along  with  the 
special  added  features  discussed  next,  enable  the  definition  of  load  paths  for 
transferring  inertial  loads  in  the  dynamics  model  to  the  structural  node  points. 
Once  these  paths  have  been  defined,  considerations  of  virtual  work  and  kinetic 
energy  lead  to  a means  for  arriving  at  rationally  determined  flexibility  and 


mass  matrices  for  the  dynamics  model.  These  latter  steps  are  fully  discussed 
in  Subsection  5. It. 


5.3. 3.1  Provision  for  Swept  Degrees  of  Freedom,  In  dynamics  modelling  it  is 
often  desirable  to  prescribe  degrees  of  freedom  that  are  parallel  to  primary 
structural  members  rather  than  the  structural  coordinate  axes.  For  example, 
pitching  motion  of  a mass  point  representing  a portion  of  the  trailing  edge 
assembly  of  a wing  might  be  described  about  an  axis  that  is  parallel  to  the 
swept  rear  beam.  To  accommodate  this  type  of  coordinate  rotation,  a feature  is 
provided  that  enables  the  program  user  to  specify  a sweep  angle,  a , at  each 
dynamic  node  point,  as  illustrated  in  Figure  5.3.  (It  is  noted  that  these  sweep 
angles  are  limited  to  the  x-y  plane. ) The  preceding  beaming  procedures  are  then 
used  to  transfer  unit  inertial  loads  in  the  local  swept  coordinate  systems  into 
structural  node  point  loads  in  the  structural  coordinate  system. 
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5. 3. 3.2  Direct  load  Transfer,  This  option  permit!  the  direct  transfer  of  force* 
and  momenta  from  the  dynamics  to  the  structure#  model  when  (elected  degrees  of 
freedom  at  a node  point  are  identical  in  both  models.  Zt  is  particularly  useful 
in  transferring  moments  to  structures  model  node  points  that  attach  beam  elements. 
Ibis  differs  frcm  the  previous  techniques  which  transfer  loads  to  node  points 
that  arc  assumed  to  be  incapable  of  sustaining  applied  moments. 

5.3.4  ftractlcsl  Considerations 

In  the  procedures  just  discussed,  applied  loadr  are  transformed  from  one  math- 
ematical model  (along  assumed  load  paths)  into  a statically  equivalent  set  of  loads 
in  the  structures  model.  Since  the  manner  in  vhich  the  user  relates  the  structural 
node  point  numbers  to  the  points  1 through  L and  N through  Q can  affect  the  final 
load  distribution,  he  should  tty  to  use  these  procedures  in  ways  that  enforce  the 
most  reasonable  load  distributions  from  the  viewpoint  of  local  structural  char- 
acteristics. This  point, is  best  illustrated  by  a simple  example. 

Consider  an  applied  wing  tip  pitching  moment  that  is  to  be  transferred  by 
four-point  beaming  to  structural  nodes  11,  12,  13,  and  lU,  and  assume  that  the 
user  has  assigned  these  node  numbers  to  I,  J,  K,  and  L as  illustrated  in  Figure 
5.4(a).  In  accordance  with  the  procedure  outlined  previously,  the  beaming  member 
e-f  will  be  essentially  horizontal,  and  the  structural  nodes  will  therefore  re- 
ceive the  applied  moment  as  approximately  vertical  forces.  Cn  the  other  hand,  if 
the  node  numbers  are  assigned  to  I,  J,  K, md  L as  illustrated  in  Figure  5.4(b), 
the  member  e-f  will  be  approximately  vertical  and  will  deliver  horizontal  forces 
to  the  structural  nodes.  If  the  structural  tip  assembly  is  attached  to  the  pri- 
mary structure  with  continuous  rows  of  fasteners  along  the  upper  and  lower  covers, 
it  would  be  more  logical  to  assign  the  node  numbers  in  accordance  with  the  second 
case. 

Some  further  recoosnendations  are  required  with  regard  to  defining  load 
beaming  from  the  dynamics  model  to  the  structures  meiel.  The  dynamics  model 
load  beaming  matrix  is  used  to  transform  the  structural  stifiPess  matrix  to 
a dynamic  flexibility  matrix  using  s procedure  discussed  in  Subsection  5.4.1. 

This  flexibility  matrix  is  subsequently  inverted  in  the  process  of  calculating 
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Figure  5-**.  Effect  of  Node  Number  Designation  on  Load  Distribution 


matrix  [B]  , which  transform  dynamics  model  displacements  to  structures  model 
displacements  (see  Equation  5.7).  The  dynamics  model  flexibility  matrix  must, 
therefore,  be  non-singular.  A singular  flexibility  matrix  can  occur  if  the 
the  number  of  structural  degrees  of  freedom  designated  for  load  beaming  is 
less  than  the  corresponding  number  of  dynamic  degrees  of  freedom  that  are 
being  created.  The  user  must  therefore  avoid  over-utilization  of  structural 
degrees  of  freedom  in  any  zone  of  the  structures  model. 

The  user  should  rlso  follow  the  general  rule  of  designating  dynamics 
nodes  which  are  adjacent  to  structures  nodes.  This  will  ensure  correct 
accounting  of  the  increment  i dynamic  mass  matrix  prescribed  by  Equation  6.2. 
Choice  of  a completely  arbitrary  dynamics  model  grid,  in  which  dynamics  loads 
are  distributed  to  a large  number  of  structures  nodes,  can  lead  to  a singular- 
ity in  the  updated  dynamic  mass  matrix  computed  by  Equation  6.1.  An  example 
of  the  recommended  procedure  is  illustrated  in  Figure  11.4  where  it  is  noted 
that  dynamic  node  points  are  coincident  in  the  X-Y  plane  with  structures  nodes 
and  are  positioned  vertically  between  the  upper  and  lower  covers. 

5.4  TRANSFORMATIONS  REQUIRED  TO  DEFINE  A DYNAMICS  MODEL 

The  following  procedure  is  * nployed  to  transform  flexibility  and  mass 
data  from  the  structures  model  the  dynamics  model.  It  assumes  that  the 
previously  discussed  beaming  procedures  have  first  been  employed  to  develop 
a transformation  matrix,  [TJ,  that  relates  forces  (or  moments)  in  the  dynam- 
ics model,  {Fjj},  to  forces  (or  moments)  in  the  structure  model  {Fs};  that  is, 

{Fs}  - [T]  {FD}.  (5.1) 

5.4.1  Flexibility  Transformation 

From  the  concept  that  virtual  work  is  invariant  under  a coordinate  trans- 
formation, it  follows  that  if  forces  relate  in  accordance  with  Equation  (5.1), 
then  displacements  relate  as 

(y  " Ct]  T K}  , (5.2) 
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vhere  f A^)  and  f Ag}  are  displacements  in  the  dynamic  tnd  structural  degrees  of 
freedom,  respectively. 

The  displacements  in  the  structure  are  related  to  applied  forces  by  the 
structural  flexibility  matrix,  that  is, 

t'Ag}  - fc.l*(F.).  (5.3) 

Substituting  (F  } from  Equation  (5.1)  into  Equation  (5.3)*  premultiplying  by 

ml 

[T]  , and  then  making  use  of  Equation  (5.2),  enables  us  to  write 

(Ap)  « Ct]  t Ck,]  -1  CtI  (fd}.  (5.4) 

This  equation  gives  dynamics  model  displacements  in  terms  of  forces  in  the 
sasie  model,  and  we  may  thus  define  the  dynamics  model  flexibility  matrix,  [a], 

as 

(A]  « [T]:PrK#]"1[T] . (5.5) 

Provision  in  made  in  FASTOP  for  the  automated  computation  of  this  flexibi- 
lity matrix  and  its  transfer  to  the  vibration  analysis  module  discussed  in  Sect- 
ion .7.  For  the  special  case,  however,  where  the  degrees  of  freedom  of  the  dynam- 
ics model  correspond  exactly  with  those  of  the  structures  model,  the  preceding 
transformation  process  may  be  bypassed.  In  tills  instance,  the  structural  stiff- 
ness matrix,  discussed  previously  in  Section  5,  is  transferred  to  and  used  di- 
rectly for  vibration  analysis. 

5.4.2  Mass  Transformation 

The  relationship  between  structural  node  displacements  and  displacements  in 
the  dynamic  degrees  of  freedom  may  be  obtained  by  substituting  Equations  (5.1) 
and  (5.4)  into  Equation  (5.3)  and  making  use  of  Equation  (5.5): 

fAg}  - [K^'Vt]  TAT^Ajj).  (5.6) 
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By  defining  the  displacement  transformation  matrix 


[BJT  - Og-HTjrAT1,  (5.T) 

and  requiring  that  the  kinetic  energy  of  the  structures  model  be  preserved  in  the 
dynamics  model,  we  may  obtain  a consistently  defined  mass  matrix  for  the  dynamics 
, model,  as  follows: 

‘ Kinetic  Energy  * l/2  (d  }^[M  (5.8) 

8 8 8 

1 where  [M  J is  a mass  matrix  corresponding  to  the  structures  model.  Substituting 

I ® 

i the  velocity  form  of  Equation  (5.6)  into  Equation  (5.8),  and  making  use  of 

| Equation  (5*7),  yields 

I * fp  m * 

i Kinetic  Energy  - l/2  l^rtEJ CMgJ [Bj  Up),  (5.9) 

( 

| where  the  inner  triple  product  is  the  desired  dynamics  model  mass  matrix,  [M^ ; 

’j  that  is, 

| 

[ [M,p  - [Bj  [Ms]  [BjT.  (5.10) 

! 

j In  the  event  that  the  degrees  of  freedom  of  the  structures  and  dynamics 

models  are  identical,  the  computation  of  [3j  may  be  bypassed,  with  the  mass 
matrix  of  the  structures  model  being  used  directly  in  subsequent  vibration 
analysis.  Further  discussion  on  the  calculation  of  the  structures  model  mass 
matrix,  in  conjunction  with  the  options  for  defining  dynamics  model  data,  is 
presented  in  Section  6. 

5.*»,3  Computational  Considerations 

Equations  (5.5)  and  (5.7)  show  that  the  inverse  of  the  structural  stiff- 
ness matrix,  [K  J”  , is  present  in  the  expressions  for  the  dynamic  flexibility 
s 

matrix  [Aj  and  the  transformation  matrix  [Bj.  As  equation  solving  routines  are 
more  eff.'cient  than  inversion  routines,  the  following  computational  procedure 
was  incorporated  into  FASXOP. 

First,  a new  matrix  [Yj  is  defined  by  the  relation 

Cy'  * t:caJ_1  [TJ,  (5.1D 

i or  equivalently, 

[KgJ  [Yj  = [Tj.  (5.12) 
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I„  vie.  .1  Equation  (5.11),  Equations  (5.5)  **  (5-7)  - »e  put  is.  the  tor. 

(4j  . [TiT  W,  <5'l3) 


and 


[Aj  [33  *=  CO 


(5.14) 


»ith,  »».  &.»>  * - “•  " irr 

(5.13)  to  compute  MS.  Then,  havins  U)  «■>!  Ms  E'>'“U°“  '5‘1 
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Section  6 


MASS  MATRIX  DEFINITION 

6.1  SUMMARY 

Two  alternate  approaches  are  available  Tor  the  computation  of  the  dynamics 
model  mass  matrix  required  for  vibration  analysis.  Ore  approach  requires  that 
initial  mass  input  data  be  specified  by  the  user  while , the  other  is  fully  auto- 
mated and  calculates  initial  mass  data  using  the  structural  model  idealization. 

6.2  INITIAL  DYNAMICS  MODEL  MASS  MATRIX  SUPPLIED  BY  THE  USER 

In  this  approach,  the  dynamics  model  mass  matrix,  [M^],  is  considered  to 
be  the  sum  of  a constant  mass  matrix,  rR^],  associated  with  the  initial  design, 
plus  a variable  mass  matrix,  which  reflects  the  accumulated  design  changes 

beyond  the  starting  design;  that  is, 

[Mp]  - ny  + ray.  (6.1) 

The  task  of  generating  [m^]  for  the  initial  design  is  usually  the  respons- 
ibility of  a weights  engineer,  who  must  account  for  the  presence  of  structural 
members,  fixed  mass  items  (equipment,  overhanging  structure,  etc.),  initial  mass- 
balance  weights,  and  all  of  the  nonoptimum  components  (fasteners,  joints,  etc.) 
which  contribute  to  the  total  weight  of  the  real  aircraft  structure.  Based  upon 
the  weights  and  locations  of  all  these  items,  the  weights  engineer  must,  develop 
a single  representative  dynamics  model  mass  matrix  for  direct  input  to  the  system. 

Computation  of  the  incremental  dynamics  model  mass  matrix,  CdMg],  is  auto- 
mated within  FASTOP.  Inasmuch  as  all  redesign  beyond  Fm^]  involves  structural 
members  and/or  mass  balances,  both  of  which  are  determined  in  the  structures 
model,  it  is  natural  to  first  compute  the  incremental  mass  matrix  in  the  struc- 
tures model,  [AM  ],  and  then  transform  the  result  to  the  dynamics  model;  that  is, 
s 

rAiy  « rs]  Tam,]'  rB]T,  (6.2) 

where  TbT  is  the  transformation  matrix  previously  defined  in  Section  5.4.2. 
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Computation  of  the  diagonal  matrix  [AM  1 ia  straightforward.  If  the  current 
gage  of  a structural  member  differs  from  its  starting  value,  the  incremental 
weight  of  that  member  (including  any  nonoptimum  factor  specified  by  the  user) 
is  computed  and  then  distributed  equally  among  its  attaching  structural  node 
points.  Bor  each  such  node,  this  "nodal  weight"  is  then  assigned  to  each  of  the 
translational  degrees  of  freedom  in  [AMg]  associated  with  that  node.  The  pro- 
cedure is  essentially  identical  for  a mass-balance  variable,  except  that  the 
incremental  weight  is  applied  entirely  to  the  single  structural  node  at  which 
the  mass  balance  is  located.  Finally,  [aMs]  ia  complete  when  all  the  structural 
members  and  mass  balances  have  been  considered. 

6.3  FULLY  AUTOMATED  COMPUTATION  OF  DYNAMICS  MODEL  MASS  MATRIX 

In  the  event  that  time  or  personnel  are  not  sufficient  to  generate  the 
initial  input  mass  data  required  by  the  preceding  approach,  a fully  automated 
mass  calculation  procedure  is  available.  At  any  stage  of  design  (starting  or 
otherwise),  the  fully  automated  method  obtains  a dynamics  model  mass  matrix, 
[Mp],  by  first  computing  a mass  matrix,  fM^],  in  the  structures  model  and  then 
transforming  the  result  by 


[Mjj]  - TB]  DTI  TB]1.  (6.3) 

The  computation  of  TM  1 is  identical  to  that  of  f AM  ] in  the  first  procedure, 
except  chat  here  total  weights  of  structural  members  and  mass  balances  ere 
used,  rather  than  incremental  weights  beyond  the  starting  design.  Also,  the 
user  has  the  option  of  inserting  fired  mass  additions  to  TM^]  (not  necessarily 
diagonal)  before  transforming  to  the  dynamics  model. 
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7.1  SUMMARY 


Section  7 

VIBRATION  ANALYSIS 


The  procedure  for  computing  modes  of  vibration  begins  with  the  transfor- 
mation of  the  familiar  structural-vibration  form  of  the  eigenvalue  problem  into 
a special  symmetric  form.  This  is  followed  by  a further  transformation  to 
tridiagonal  form  by  Householder's  method  (Reference  7-l>  page  290).  A Sturm 
sequence  technique  (Reference  7-1,  page  300)  is  then  employed  to  determine  the 
eigenvalues,  and  inverse  iteration  (Reference  7-1,  page  321)  is  used  to  cal- 
culate the  eigenvectors.  This  procedure  has  proved  to  be  both  efficient  and 
accurate  for  problems  where  the  solution  may  be  achieved  directly  in  core. 

7.2  TRANSFORMATION  OF  EIGENVALUE  PROBLEM  TO  A SPECIAL  SYMMETRIC  FORM 

The  first  step  in  the  solution  procedure  is  to  transform  the  structural- 
vibration  eigenvalue  problem  into  the  form 

[D]  fy}  * XfY},  (7.1) 


where 


[D 1 = a real  symmetric  matrix  having  only  real  roots 

fY)  * an  eigenvector  of  the  transformed  problem 

X a a real  eigenvalue  (=  l/oi  ) 

id  = a natural  frequency  in  radians/sec. 

This  transformation  can  take  either  of  two  forms  depending  on  whether  the 
initial  structural  representation  is  in  terms  of  a stiffness  matrix  or  a 
flexibility  matrix.  It  should  be  noted  that  the  eigenvalue,  X,  is  defined 
here  as  l/<£  since  the  solution  procedure  to  be  defined  subsequently  will 
determine  the  higher  eigenvalues  with  greatest  accuracy. 

7.2.1  Stiffness  Matrix  Formulation 

The  eigenvalue  problem  for  the  stiffness  matrix  formulation  may  '->e  . 
written  as 


^5 


Cm]  fx]  = xfKl  [x], 


(7.2) 


where 

fM]  * the  mass  matrix  of  the  dynamics  model  idealization 

[K]  = the  stiffness  matrix  of  the  dynamics  model  idealization 

fx)  » an  eigenvector  in  dynamics  model  coordinates 

and  the  definition  of  X is  the  same  as  for  Equation  (7.1).  Factoring  the 
stiffness  matrix  into 


[K]  - [L]  fL)T  (7.3) 

by  Cholesky  decomposition,  discussed  previously  in  Section  4.3.1,  and  sub- 
stituting into  Equation  (7.2)  yields 

CM]  fx]  = X[L]  [if  fx],  (7.4) 

where  the  factorization  of  Equation  (7*3)  requires  that  [K"l  be  nonsingular. 
Letting 


{y}  = [L]T  fx]  (7.5a) 

or 

Cx]  = [L]”T  fY},  (7.5b) 

and  premultiplying  Equation  (7.4)  by  [L]"^  gives 

[L]"1  Cm3  [l]-t  Cy]  = x fy)  (7.6) 


or 


Id]  Cy)  = x (y},  (7.7) 

where  [d]  = [l]"^  [m]  Cl>]”^.  It  should  be  noted  that  [d]  is  symmetric  since 
Cm]  is  always  symmetric. 
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Alter  the  eigenvalue  problem  of  Equation  (7.?)  is  solved,  as  discussed 
subsequently,  the  eigenvectors  are  transformed  back  to  the  dynamics  model  co- 
ordinate system  using  Equation  (7.5b). 

7.2.2  Flexibility  Matrix  Formulation 

The  eigenvalue  problem  for  the  flexibility  matrix  formulation  is 

[Ml  fx)  * X (A]'1  fxJ  (7.8) 

or 

[A]  [M]  fx3  = X fx},  (7.9) 

where  [Al  = [K]'1  * the  flexibility  matrix,  and  Tm3,  [K],  fx),  and  X are  the 
same  as  defined  previously.  Factoring  the  mass  matrix  into 

[M]  = [L]  [L]T  (7.10) 

and  substituting  into  Equation  (7*9)  yields 

TA]  [L]  -Lf  fx}  = X fx}.  (7.11) 

Making  the  same  substitution  for  fx}  as  defined  by  Equation  (7.5b),  but,  in  this 

T 

case  using  [L]  as  defined  above,  and  premultiplying  by  [Ll  , yields 

[l]t  [a]  [Li  fy)  - x fy}  (7.12) 

or 

[Dl  fy}  = x fy},  (7.13) 


where 


TdI  - [bf  [A]  [L). 

Here  again,  symmetry  is  preserved  in  the  transformation,  and  the  eigen- 
vectors are  transformed  by  Equation  (7.5b). 
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7.3  COMPUTATION  OF  VIBRATION  MODES  FOR  A CANTILEVER  STRUCTURE 

Computation  of  vibration  modes  for  a cantilever  structure  proceeds  as 
described  previously  in  Section  7,  where  the  formulation  of  the  eigenvalue 
problem  employs  either  the  structural  stiffness  matrix  (Equation  7-2)  or  the 
dynamics  model  flexibility  matrix  (Equation  7-9)  of  the  supported  structure. 
The  former  approach  is  taken  where  degrees  of  freedom  of  the  dynamics  model 
correspond  exactly  with  those  of  the  structures  model. 

7.4  COMPUTATION  OF  VIBRATION  MODES  FOR  A FREE-FREE  STRUCTURE 


Although  the  computation  of  free-free  modes  also  uses  the  stiffness  or 
flexibility  matrix  of  the  supported  structure,  the  formulation  of  the  free- free 
eigenvalue  problem  requires  relaxation  of  the  fixed  support  points  to  allow 
rigid-body  motion  of  the  structure.  The  analytical  procedures  are  described 
below.  In  the  following  discussion  a hypothetical  "plug"  is  defined  which  is 
assumed  to  be  rigidly  interconnected  to  all  the  fixed  support  points  of  the 
structures  model;  the  plug  mass  properties  represent  that  portion  of  the  config- 
uration not  included  in  the  dynamics  model  of  the  flexible  structure. 

Figure  J.l  shows  an  unsupported  (free-free)  configuration  consist* ng  of 
a rigid  "plug"  section  connected  at  points  A and  B to  flexible  structure; 
local  flexibilities  can  exist  at  A and  B. 


Typical  Dynamic  Point  on  Flexible  Structure 


The  vibration  equation  for  this  unsupported  system  may  be  written  as: 


} 


(7.14) 


where 


M “ a vector  of  absolute  displacements  of  dynamic  points  on  the 
flexible . structure . 


is  a vector  of  absolute  displacements  of  the  plug  reference  point. 

is  the  dynamic  mass  matrix  of  the  flexible  structure  alone;  this 
matrix  was  discussed  in  Sections  6.2  and  6.3. 

is  the  mass  matrix  of  the  plug  alone. 

is  the  stiffness  matrix  for  dynamic  degrees  of  freedom  of  the 
flexible  structure  excluding  the  plug;  that  is,  this  matrix  defines 
forces  at  dynamic  points  on  the  flexible  structure  due  to  displace- 
ments at  those  points  alone. 


is  the  stiffhess  matrix  for  the  plug  degrees  of  freedom;  that  is, 
this  matrix  defines  the  forces  acting  on  the  plug  due  to  motion  of 
the  plug  alone. 


M 

[*»] 


is  the  stiffhess  matrix  defining  forces  at  dynamic  points  on  the 
flexible  structure  due  to  plug  motion  alone. 

in  the  stiffhess  matrix  defining  forces  on  the  plug  due  to  motion 
of  the  flexible  structure  alone. 


Now,  the  absolute  motion  of  the  flexible  structure  {^d}  is  due  to  both  rigid- 
body  motion  of  the  entire  configuration  and  flexible  (or  relative)  motion, 
I^DR  } ’ Accordin«1y>  if  the  motions  are  chosen  to  prescribe  the  rigid- 
body  motions,  it  follows  that 


(7.15) 
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where  [^DJ  defines  the  displacements  of  the  dynamic  points  on  the  structure 
due  to  unit  rigid-body  motions.  The  basic  vibration  equation  can  now  be  recast 
in  terms  of  i^DR^  and  l^p}  by  first  substituting  Equation  (7*15)  into  Equation 
(7- 1^0  and  then  premultiplying  the  result  by  the  transpose  cf  the  transformation 
matrix  in  Equation  (7.15).  The  result  is 


Tad  iTTj  ) (Tad  i a J( 


where 

W - W U . W ■ kf  h] 

w-krwy-w 


and 

M x [*0]  [xd]  + [^p] 

n = yT  w 


(7.16) 


(7.17a) 


(7.17b) 


However 


, [kd»]  , [kad]  «a  H 


must  all  be  zero  because  rigid-body  motion  alone 


(as  prescribed  by  {^p}  ) cannot  induce  restoring  forces  from  the  stiffness 
matrix  of  the  free-free  system.  The  same  conclusion  can  be  reached  in  another 
manner.  Consider,  for  example,  the  two  terms  present  in  the  expression  for 


w 


in  Equation  (7.17b).  These  two  terms  respectively  define  the  forces  acting 
at  dynamic  points  on  the  structure  due  to  (a)  rigid- body  motion  of  the  flexible 
structure  alone  (see  Figure  7«2a),  and  (b)  motion  of  the  plug  alone  ^see  Figure 
7.2b).  As  these  two  force  systems  are  negatives  of  each  other 


, |>A 


must 


vanish.  Matrices  [KAdJ  and  J^aJ  can  be  shown  to  vanish  in  a similar  manner. 
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* UD] 


Figure  7*2.  Rigid-Body  Motions  of  Structure  Alone  and  Plug  Alone 


Thus,  Equation  (7.16)  may  be  rewritten  as 


U) 


k_,MJ  + Sj.® 

pAD  J*]  ppj  L5]0] 


*0 


(7.18) 


It  follows  from  Equation  (7.18)  that  the  rigid-body  coordinates  { } are 
expressable  in  terms  of  the  relative  coordinates  { ^DP.f  ; that  is 


l^p  I - - [ma  ] [mad]  {^dr}  . 


(7.19) 


Finally,  if  Equation  (7.19)  Is  substituted  into  Equation  (7.18),  the  vibration 
equation  for  the  free-free  structure  in  terms  of  relative  coordinates  becomes 

- u,2  [MT)Ff]  {%<}  + »{o}  , (7*20) 

where 

[MDFf]  * ft]  " [%]  [Ma]  [MAd]  . (7.21) 
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Equation  (7.21)  defines  the  n&ss  matrix 


^dffJ 


for  a free-free  vibration 


enalysis  regardless  of  vhieh  of  the  two  alternate  approaches  is  used  to  obtain 
the  dynamics  model  mass  matrix  M (see  Sections  6.2  and  6.3).  Of  course, 
the  plug  mass  matrix [H*!  must  be  supplied.  If  the  flexibility  approach  is  U3ed, 


Equation  (7.20)  may  be  written  in  the  form 


p[>]  [>]  >['])  j t ■ “•  <7  S2: 

where  is  the  flexibility  matrix  for  dynamic  degrees  of  freedom  of  the 

flexible  structure  excluding  the  plug  (see  Equations  7.8  and  7. 9)* 

After  the  {^Dh|  are  generated  in  the  solution  of  Equation  (7-20)  or 
Equation  (7.22)  using  the  eigenvalue  solution  procedures  described  previously, 
the  associated  plug  motions  can  be  obtained  from  Equation  (7.19)  and  the 
absolute  displacements  of  the  dynamic  points  on  the  structure  may  be  confuted 
from  Equation  (7.15).  The  generalized  mass  matrix  for  the  free-free  system 
may  then  be  obtained  from  the  relation 

[ 

In  the  event  that  the  degrees  of  freedom  of  the  structures  and  dynamics 
model  are  identical,  matrix  [*d]  is  replaced  by  matrix  D's]  which  defines  the 
displacements  of  the  structures  points  due  to  unit  rigid-body  motions.  The 
matrix  [xj  is  always  required  if  flutter-velocity  derivatives  are  to  be 
computed  for  a free-free  vibration  model  (see  Subsection  9*3.1). 


«] 


0 j" pj  \0p j 


(7.23) 


Section  8 


FLUTTER  ANALYSIS 


8.1  SUMMARY 

This  module  of  the  FASTOF  package  determines  the  oscillatory  pressures 
and  the  generalized  aerodynamic  forces  for  the  lifting  surface  to  be 
analyzed,  given  a set  of  normal  mode  shapes  and  frequencies.  Additionally, 
given  the  generalized  masses  corresponding  to  the  modes,  the  program  solves 
the  flutter  equation  to  determine  the  flutter  speed  and  values  of  modal  damp- 
ing  and  frequency  as  functions  of  air  speed. 

Generalized  aerodynamic  forces  are  computed  using  either  the  subsonic 
assumed-pressure-function  procedure  (kernel  function),  tne  supersonic  Mach- 
box  method  or  the  subsonic  doublet-lattice  procedure.  In  determining  the 
flutter  speed,  these  aerodynamic  forces  are  required  at  many  different  re- 
duced frequencies;  consequently,  to  save  computing  time,  these  forces  are  j 

determined  at  the  required  reduced  frequencies  by  interpolation,  using  a <’ 

small  number  of  directly  calculated  aerodynamic  forces  as  a basis.  \ 

The  flutter  solutions  are  obtained  by  use  of  either  the  conventional  fc. 
method  or  an  improved  version  of  the  p-k  method  of  Reference  8-1.  To  allow 
the  user  to  study  the  flutter  mechanism,  several  parameter  variations  have 
been  automated.  For  redesign  purposes,  the  eigenvectors  and  their  associated  \ 

row  vectors  are  determined  at  the  flutter  speed  when  the  p-k  method  is  \ 

selected. 


/ 

/ 


/ 
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8.2  AERODYNAMICS  ROUTINES 

Three  oscillatory  lifting-surface  aerodynamics  routines,  each  based  on 
linearized  potential  flow  theory,  are  included: 

(1)  a subsonic  assumed-pressure-function  program  (kernel  function)  for 
analyzing  one  or  more  planar,  noninteracting,  lifting  surfaces; 

(2)  a supersonic  Mach-box  program  for  analyzing  one  or  several  planar, 
noninteracting,  lifting  surfaces; 

(3)  a subsonic  doublet -lattice  program  for  analyzing  a general  con- 
figuration of  multiple  interacting  nonplanar  surfaces  and  slender 
bodies. 

To  minimize  the  required  input  data  preparation,  all  of  these  programs  have 
been  written  or  modified  to  provide  automated  geometry  definition  and  modal 
interpolation  of  vibration  data.  The  latter  either  may  be  obtained  on  a 
magnetic  tape  (or  disk)  from  the  vibration  analysis  program  described  in 
Section  7 or  may  be  user-supplied  on  input  data  cards.  With  the  former 
method  of  providing  the  vibration  data,  the  user  has  the  option  of  elimi- 
nating modal  data  not  needed  in  the  flutter  analysis.  To  conserve  computer 
machine  time,  the  program  can  be  instructed  to  save  the  aerodynamic  influence 
coefficients  on  magnetic  tape  for  subsequent  reanalysis  with  altered  vibra- 
tion modes. 

The  following  is  a brief  theoretical  discussion  of  the  aerodynamics 
routines  and  the  associated  options. 

8.2.1  Subsonic  Assumed-Pressure-Function  Program 

In  the  subsonic  regime,  aerodynamic  forces  are  computed  using  the  assumed- 
pressure-function  method  of  Reference  8-2.  For  the  general  nonplanar  case, 
the  pressure  on  a harmonically  oscillating  surface,  such  as  shown  in  Figure 
8.1,  is  related  to  the  downwash  by  the  following  integral  equation  derived 


I 

i 

i 


U is  the  free  stream  velocity 

h,(*  are  the  normal  displacement  and  the  streamwise  slope  of  the  surface 
ui  is  the  harmonic  frequency  of  oscillation 
p is  the  air  density 
t is  time 

K is  the  kernel  function  specifying  the  normalvash  angle  at  P^  due 

to  a unit  pressure  at  P 

k is  the  reduced  frequency  * uib  /U 

b is  the  reference  semichord 
o 

M is  the  Mach  number. 

This  can  be  rewritten  as 


"(P,)  * j-~-  Jf  p(P)K(Pj,  P,  k,  M)d§do 


"(Pj)  * ^Jfcp(V)K(?y  P,  K,  M)d$da, 


(8.2) 


and  for  planar  surfaces  can  be  written  as 


»(x,y)  ’ 8?  ffcp(g,Tl)K(x-5,  y-D,k,M)d?dTl, 


where 


w is  the  normal  or  downwash  angle  a w/u 

Cp  is  the  differential  pressure  coefficient  * p/q 

q is  the  dynamic  pressure  * ^ pU2. 

In  the  program,  the  planfora  coordinates  (x,y)  and(5,T|)  are  normalized 
with  respect  to  the  root  semichord,  bc.  These  nondimens ional  coordinates 
(x,y)  and  are,  in  turn,  transformed  into  the  coordinates  of  a square 

planform,  (x,  y)  and  (§,T))  as  shown  in  Figure  8-2.  After  transformation 
of  coordinates,  the  integral  equation  becomes 

w(x,y)  « — + ikh) 


, — 2_  f r c (?,Tl)K(x-t,  y-l[,k,M)b(Tl)dCdTl, 

8*o  JJ  


Figure  8.1.  Nonplanar  Harmonically  Oscillating  Surface, 
from  linear  potential  theory: 

*(Pj)  *eiUlt  * 5^-  jfv^Wy  P,  M)  d^da, 

s (8.1) 

where 

P is  any  point  on  the  surface,  the  coordinates  of  which  are  §,  T),  Q 
P.  is  a point,  the  coordinates  of  which  are  x,  y,  z 

x, §  are  streamwise  coordinates 

y, Tl  are  spanwise  coordinates 

z, C  are  vertical  coordinates 

a is  the  tangential  spanwise  coordinate 

(equivalent  to  T|  for  a planar  surface) 

S is  the  total  lifting  surfaje  area 

p is  the  pressure  difference  between  the  upper  and  lower  covers  of  the 
surface 

w is  the  complex  normal  wash  = U a + iuh 


t 


U is  the  free  stream  velocity 

h,a  are  the  normal  displacement  and  the  streamwise  slope  of  the  surface 
o)  is  the  harmonic  frequency  of  oscillation 
p is  the  air  density 
t is  time 

K is  the  kernel  furction  specifying  the  normalwash  angle  at  Pj  due 

to  a unit  pressure  at  P 

k is  the  reduced  frequency  * ubo/U 

b is  the  reference  semichord 
o 

M is  the  Mach  number. 

This  can  be  rewritten  as 

"(Pj)  « P>  k’  M)d5d0 

s 


v(pj)  * (fcp^K(py  p»  K)  M)d?do»  (8-2) 

s 

and  for  planar  surfaces  can  bo  written  as 

v(x,y)  « gi-  jjJcp(5,Tl)K(x-5,  y-T1,k,M)d?dTl,  (8.3) 

S 

where 

w is  the  normal  or  downwash  angle  *>  w/U 

Cp  is  the  differential  pressure  coefficient  * p/q 

q is  the  dynamic  pressure  * b pU2. 

In  the  program,  the  planform  coordinates  (xyy)  ar.d  (5 » "H)  are  normalized 
with  respect  to  the  root  semichord,  b0.  These  nondimensional  coordinates 
(x,y)  and  (^,“71)  are,  in  turn,  transformed  into  the  coordinates  of  a square 
planform,  (x,  y)  and  (5,T))  as  shown  in  Figure  8-2.  After  transformation 
of  coordinates,  the  integral  equation  becomes 

"(x,y)  = C-  + ikh) 

o x 


-1  -1 
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Planform  (in  real  coordinates) 


Transformations  to  Nondimensional  Coordinates 
x = x/bQ  5 > 5/b0 

7 - yA>0  * = V*>0 

Transformations  to  Normalized  Nondimensional  Coordinates 

_ x-xm(y)  - = 5-5m(Tl) 

- ~ b(y)  - b('Ti)  ' 

y - y/i0  2 » T)/10 

The  above  transformations  convert  the  planform  into  a unit  square. 
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where 


1 is  the  semispan 
o 

b is  the  root  semichord 
o 

b(T))  is  the  local  semichord. 

In  the  above  equation,  the  planform  has  been  assumed  symmetrical  about  y * o 
as  indicated  by  the  limits  of  spanwise  integration.  Equation  (8.4)  must 
now  be  solved  for  the  unknown  pressure  coefficient  distribution,  Cp  (5,3) , 
in  terms  of  the  known  dcrwnwash  distribution.  The  solution  procedure  requires 
that  the  pressure  distribution  be  expressed  as  a series  of  terms,  each  of 
which  is  the  product  of  an  unknown  constant  coefficient  and  a specific  load- 
ing function.  The  unknown  constant  coefficients  are  then  moved  outside  the 
integral  and  the  remaining  integration  is  performed.  To  assure  that  the 
series  can  represent  the  pressure  correctly,  the  loading  functions  are 
chosen  to  satisfy  the  following  boundary  conditions: 

• A square-root  singularity  at  the  planform  leading  edge 

• A zero  at  the  trailing  edge 

• A zero  at  the  planform  tip  (or  outboard  edge) 

• An  infinite  slope  at  the  tip. 

A suitable  3eries  expression  satisfying  these  conditions,  which  is  used  in 
Reference  8-2  and  in  this  program,  is 

CP  (§’3>  * STn) 

1 

I . 

•Ja1  J ^ J ai  1 - (l-+§)  I > 

V1*3  ri  1,3  ■ / (8.5) 


where 


a.  . are  the,  as  yet,  undetermined  pressure  series  coefficients 

0 = 0 for  a symmetrical  distribution  of  airloads  about  the  root  and 
= 1 for  an  antisymmetrical  distribution 

1 is  the  number  of  terms  in  the  chordwise  series 
J is  the  number  of  terms  in  the  spanwise  series. 


After  substituting  the  series  expression  for  the  differential  pressure 
into  Equation  (8.4)  and  removing  the  unknown  constants  from  the  integral, 
the  remaining  chordwise  integrations  are  of  the  form: 

1 

Id  « J*  (l-u)ar(l+u)Sfd(u)du,  (8.6) 


where 

u “ 5 

a * ^ and  8 » 

fd(u)  is  an  appropriate  polynomial  of  degree  d. 

In  the  spanwise  direction,  the  integrals  also  are  of  the  above  form  but  with 
different  values  of  a and  8 . Inspection  of  Equation  (8.5)  shows  the  lead- 
ing spanwise  factor  to  be  Vl  - if;  however,  the  kernel  function,  K,  of 
Equation  (8.4)  contains  a factor  of  l/(l  - T^).  Consequently,  the  product 
of  the  pressure  series  and  the  kernel  function  has  a factor  of  l/^l  - Tp 
making  the  appropriate  values  of  a and  8 *•  of  = 8 =»  in  the  spanwise 
integral. 

These  integrals  are  evaluated  using  the  Gauss-Mehler  quadrature  of 
Reference  8-3  (pages  312-357)  to  give 

id-i:  "KKV*  (8-7) 

q=l 

where 

Q,  is  the  number  of  integration  points 

H are  weighting  functions  given  in  rows  two  and  four  of  Table  8.1 
u^  are  integration  points  given  in  rows  two  and  four  of  Table  8.1. 

With  this  method,  if  the  integrand  can  be  exactly  represented  by  the  product 
of  a simple  singular  function,  such  as  and  a polynomial  of 

degree  d = 2Q-1,  then  using  the  Q integration  points  would  produce  the 
integral  with  no  error.  Although  the  true  integrand  cannot  be  represented 
exactly  with  a finite  number  of  terms  in  the  polynomial  part,  by  choosing 
enough  integration  points,  the  error  may  be  made  S3  small  as  desired. 
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TABLE  8.1.  GAUSS -MEHLER  QUADRATURE  FORMULAE 
1 Q 

For  j (l-u)tt(l+u)Pfd(u)du  » ^H(uq)fd(uq) 

-1  q-1 


or 

8 

H ( weighting  function) 

Uq  (integration  point) 

£ 

£ 

m <lV> 

* cos  (?5i> 

£ 

-£ 

2Q+1  ^”Uq^ 

- 008  (iqTr  n) 

-£ 

£ 

- COS  (gL) 

-£ 

-£ 

TT 

5 

- C08  n) 

After  the  Integrations  are  performed,  Equation  (8.4)  can  be  written 

v(Z,y)-£  Y,  Li,J(?’7)  *i,j ' 

i*i  j«i 
N 

■y;  Ln(x,yO  &n,  (8.8) 

n=l 


where 


*i>J’  N 


IxJ 


the  result  of  the  integration  of  the  product  of  the 
kernel  function  and  the  i,j^  pressure  series  term. 


When  Equation  (8.8)  is  applied  at  a set  of  points  on  the  wing  (called 
collocation  points),  a resulting  system  of  linear  equations  relating  the 
downwash  at  these  points  to  the  unknown  coefficients  is  formed: 

N 

w(P  ) = >L  (P  ) • a 
c'  nv  c'  n 

n=l 


) 


J 


s 
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> 
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t 


s [L]  M 5 (8.9) 

(CxN)  (Nxl) 

where 

Pc  is  the  c 1 collocation  point 
C is  the  total  number  of  collocation  points. 

The  locations  of  the  collocation  points  are  chosen  so  that  the  error 
that  would  result  from  calculating  the  lift  based  on  fitting  a polynomial 
thru  the  collocation  points  is  minimized.  Reference  8-2  shows  that  in  the 
chordwise  direction,  this  loading  is  calculated  from 

1 

f(x)  dx,  (8.10) 


or 


\v\ 

(Cxi) 


where 
of  the 


f(x) 

form 


is  a polynomial.  In  the  spanwise  direction,  the  integral  is 


g(y)  dy, 


(8.11) 


where  g(y)  is  a polynomial.  Since  these  integrals  are  of  the  same  form  as 
Equation  (8.6),  Gaussian  quadrature  can  again  be  used.  In  the  streamwise 
case.  Equation  (8.10),  a = and  0 = +|,  while  in  the  spanwise  case  of 
Equation  (8.1l)  cr  * 8 » +J.  Hence,  the  proper  choice  of  collocation  points 
can  be  found  in  the  first  and  third  rows  of  Table  8.1. 

The  kernel  function,  K,  of  Equation  (8.4)  possesses  a singularity  as 
y •*-tT  and  x'-*'^.  Consequently,  the  collocation  and  integration  points  must 
not  be  allowed  to  coincide.  From  the  above  discussion  and  Table  8.1,  the 
chordwise  locations  of  these  points  are: 

pj  i 

integration  points  vl^  = -cos  , i = l,...,Ix 

collocation  points  uc  = -cos  (|g-“1)  c * l,...,Cx,  (8.12) 

where 

lx  is  the  total  number  of  chordwise  integration  points 
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is  the  total  number  of  chordvise  collocation  points. 
In  the  spanwise  direction,  the  locations  are: 

2i-l 

integration  points  * -?os  (^~—  n) , i * l,...,Iy 

y 


collocation  points  u 


-cos  <§=*>,  C 


(8.13) 


wnere 


I is  the  total  number  of  spanwise  integration  points 

y 

C is  the  total  number  of  spanwise  collocation  points. 

With  a choice  of  Iv  » C +1  as  in  Reference  8-2,  the  points  are  distinct.  In 

y y 

the  current  program,  ly  cc..  be  chosen  differently  as  -will  be  described  later. 

Equation  (8.9)  can  be  written  for  each  mode  to  be  used  in  the  analysis; 
hence,  for  several  modes 

[w]  » £aj  , (8.l4) 

(CxNM)  (CxN)  (NxNM) 

where  NM  is  the  total  number  of  modes.  In  Reference  8-2,  the  total  number 
of  collocation  points,  C,  must  equal  the  total  number  of  undetermined  pres- 
sure series  coefficients,  N.  With  such  a restriction,  is  a square 
matrix  and  Equation  (8.l4)  can  be  solved  for  [a]  by  simple  inversion.  For 
relatively  complicated  mode  shapes,  requiring  a large  number  of  collocation 
points,  this  restriction  forces  the  use  of  high  order  polynomials  to  repre- 
sent the  pressure  distribution.  Since  this  can  lead  to  unrealistic  convolu- 
tions in  the  calculated  pressures,  a new  procedure  has  been  adopted  in  the 
present  program  which  permits  the  use  of  fewer  polynomial  terms.  In  this 
approach,  the  [l]  matrix  is  no  longer  square  ^ C>N  ^ and  Equation  (8.l4) 
is  solved  in  a least-square  sense  as  described  in  Reference  8-4. 

A least-square  solution  to  Equation  (8.l4)  consists  of  a matrix  [a] 
which  minimizes 


[w]  - fL][a] 


(8.15) 


where 


indicates  the  Euclidean  norm  - defined  rs  the  square  root  of 
the  sum  of  the  squares  of  the  terms  of  the  array.  This  is  equivalent  to 
minimizing  ( [„]  I _ ^ I (l)  . [t]  [ .]  ).  (8.l6) 

for  each  of  the  NM  columns  of  [w]  - [L][a]. 


62 


When  the  rank  of  J equals  N,  a solution  can  be  obtained  hypothetically  by 


"1 


m 


(8.17) 


H[*Hf 

■{['my 

However,  (lJ  T [l]  is  frequently  ill-conditioned  making  this  direct  approach 
impractical.  Reference  8-4  shows  that  thi3  problem  can  be  avoided  by  first 
decomposing  [l]  by  an  orthogonal  transformation  [t]  , formed  by  the  continued 
product  of  Householder  transformations  (Reference  8-5),  such  that 

HM  -[£]• 

(CxC)(CxN)  (CxN) 


(8.18) 


where  [tfj  is  an  Nxll  upper  triangular  matrix.  Applying  this  transformation 
to  both  sides  of  Equation  (8.l4),  one  obtains 


[-!-][•] -W» 

Upon  extraction  of  the  first  N rows,  Equation  (8.19)  becomes 

[r]  [ a ] = [l'0][T][wj, 


(8.19) 


(8.20) 

U J 

where  CO  is  an  NxN  identity  matrix.  Since  R is  an  upper  triangular  matrix, 
this  equation  is  easily  solved  by  the  back  solution  part  of  any  linear- 
system-solver  (see  Reference  8-3,  pp  428-429).  It  can  be  shown  that 
the  solution  of  this  equation  is,  indeed,  a least  squares  solution  to 
Equation  8.14.  This  is  done  by  showing  that  it  can  be  reduced  to 
Equation  8.17.  To  reduce  the  computing  time  required  in  subsequent 
re-analyses  made  with  revised  modal  data,  the  user  can  save  the  matrices 
[tJ  and  [t| • [l]  for  future  use.  In  subsequent  analyses,  new  [w] 
matrices  are  generated  and  a system  analogous  to  Equation  8.20  is 
formed: 


[Ho]  * ([t]  • [l])  .fa]  = [l ; 0 1 - [ T ] • [w]  , (8.21) 

L J eld  new  old  new 


from  which  pressure  series  coefficients,  [a],  can  be  determined  for  the 
new  modes. 


1 


•1 
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Once  the  pressure  series  coefficients  have  been  determined,  the 
normalized  generalized  air  forces  are  computed  by 


— ft 
2 Hrs 


(x,y)  P8(x,y)  dx  dy, 


(8.22) 


where 

ft  is  the  generalized  air  force  in  the  r mode  due  to  the  pressure 
is  th 

from  the  s mode 

S/2  indicates  integration  over  one  half  of  the  planform 

h is  the  deflection  distribution  in  the  r^  mode 
r iu 

p is  the  pressure  distribution  in  the  s mode. 

s 

This  can  be  rewritten  as 


(x,y)  C (x,y)  dx  dy. 
ps 


(8.23) 


Transforming  to  nondimensional  coordinates  and  substituting  Equation  8.5, 
this  becomes 


rs 


dsli  £ jf? 


„2(j-l)  +« 


-1 


a^j^x+ljx1-2^  h^(x,y)dxdy. 

(8.2b) 


This  integration  is  performed  using  the  appropriate  Gauss-Mehler  quadrature 
formulae  previously  described.  The  number  of  integration  points  used  in  the 
chordwise  integration  is  twice  the  number  of  chordwise  collocation  points,  c , 
and  in  the  spanwise  integration  equals  the  number  of  collocation  points,  C^. 


Important  considerations  in  using  the  assumed-pressure  function  program 
are  the  choices  of  the  number  of  terms  in  the  pressure  polynomial,  the  number 
of  integration  points,  and  the  number  of  collocation  points.  For  Equation  8.1b 
to  have  a unique  solution,  the  number  of  collocation  points  in  both  the 
chordwise  and  spanwise  directions  must  be  equal  to  or  greater  than  the  num- 
ber of  polynomial  terms  chosen  to  represent  the  pressure  distribution  in 


6b 


the  respective  directions;  i.e., 


Cx  * I and 

Cy  » J.  (8.25) 


Furthermore,  for  chordwise  and  spanwise  integrations  of  the  product  of  the 
pressure  and  the  kernel  function  to  be  satisfactory,  Reference  8-2  advises 
that  the  number  of  integration  points  in  each  direction  should  obey  the 
relationships 


Ix  s Cx  and 


ly  * Cy+1. 


(8.26) 


Tc  increase  accuracy,  more  integration  points  can  be  used  in  the  present 
program  than  the  minimum  recommended.  In  Reference  8-6,  it  is  shown  that 
the  excess  chordwise  integration  points  can  cause  numerical  difficulties 
if  their  number  is  not  chosen  carefully  to  avoid  close  proximity  between  the 
collocation  and  integration  points.  It  suggests  the  following  formula  to 
govern  the  choice: 


Ix  - (Cx  + i)  • .(2NC-1), 


(8.27) 


where  NC  is  a positive  integer  and  I is  taken  as  the  truncated  integer 
value  of  the  expression.  In  the  present  program,  the  user  selects  and  NC 
as  data  after  which  the  program  itself  computes  from  the  above  equation. 
Additionally,  the  user  supplies  Cy  and  a positive  integb.  , NS,  as  data;  and 
the  program  computes  the  number  of  spanwise  integration  stations  by  the 
formula: 


Iy*NS(Cy+l).  (8.28) 

One  final  empirical  guideline  suggested  in  Reference  8-6  to  obtain 
converged  results  is  to  choose  the  number  of  collocation  points  to  satisfy 
the  ratio 


Cx  ~ V'T-M2  (A.R.) 
<y  ~ cos  An 


(8.29) 


where 

M is  the  Mach  number 
A.R.  is  the  aspect  ratio 

Affl  is  the  sweep  of  the  wing  midchord  at  the  root. 


where 


x,  § ere  streamwise  coordinates 

y,  T|  are  spanwise  coordinates 

p is  the  differential  pressure  between  the  upper  and  lower  covers 
of  the  surface 
p is  the  air  density 

U is  the  free  stream  velocity 

u>  is  the  frequency  of  oscillation 

<p  is  the  velocity  potential 

S is  the  lifting  surface  area  bounded  by  the  inverse  Mach  cone 
emanating  from  (x,  y) 

w is  the  complex  downwash  velocity  s U a + ituh  * Uw 
h,  a are  the  deformation  and  slope  of  the  surface 
M is  the  Mach  number 

a * y^-i  

r « V(*-§)2  " p2  (y-1))2 

is  the  differential  pressure  coefficient  * p/J  pU2. 

With  the  exception  of  special  cases,  the  integral  cannot  be  evaluated  in 
closed  fora;  hence,  a numerical  approach  is  required.  In  Reference  8-7,  the 
area,  S,  is  divided  into  elementary  small  rectangular  boxes  having  their 
diagonals  parallel  to  the  Mach  lines  as  shown  in  Figure  8.3-  The  rectangles 
are  subsequently  converted  to  squares  through  the  coordinate  transformations 

Streamwlse  ST  * x/b, 

? = 5/b, 

Spanwlse  y = £y/b 

• “H  = eVb,  (8.31) 

where  b is  the  streamwise  dimension  of  a box  and  b/3  is  the  spanwise  dimension. 
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Assuming  the  downwash  is  constant  ovc  each  of  these  "Mach  boxes,"  Equation 
(8.30)  can  he  rewritten  as: 


Cp  (x,  y ) 


•»]•//•■'«' -t 


k R 


.2 


M b ~ 

— — d§  d7)  , 

b p 


s&»-k2:v&  *“>//•■’ 

d 


ikM  COS^.liB, 

M R 


3 


or  Cp  (x,  y)  - Yj  v3  Cj(*’  30, 

d 


(8.32) 

(8. 33) 


where 


S3  is 

..  .th 
tne  3 

box  within  the 

inverse  Mach  cone  emanating  from  (x,  y) 

k = 

bu 

U 

»v 

k = 

(»)% 

*s* 

R » 

{(T- 

?)  - (y  -TU 

b 

is  the 

streamwise  box 

size 

Wj 

is  the 

downwash  on  the  j box 

C is  the  3 pressure  influence  coefficient  for  point  (x-  y),  i.e. , 

J th 

the  pressure  at  the  point  due  to  a unit  downwash  on  the  J box. 


Two  metitods  have  been  classically  used  to  perform  the  complex  integration 
over  the  box  area.  The  first,  developed  in  Reference  8-8,  uses  a mean  value 
of  the  exponential  and  cosine  terms  when  a box  is  far  removed  fran  the  point 

.is/  'w.  2 

(x,  y)  >d  a series  expansion  up  to  k when  a box  is  close.  Although  these 
approximation:-  simplify  the  integration,  they  introduce  significant  errors, 
when  the  reduced  frequency,  k,  is  high.  For  a more  exact  evaluation,  a second 
method,  a Bessel  function  series  representation  of  the  integral,  i.e  presented 
in  Reference  8-7. 
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The  method  used  in  the  current  program  is  different  from  both  of  these 
two  in  that  Gaussian  quadrature  is  used  to  evaluate  the  integral.  With  this 
technique,  singularities  that  occur  in  the  integrand  when  the  box  area  is  cut 
or  touched  by  the  inverse  Mach  cone  emanating  from  the  point  (x , y)  - see 
Figure  8.4a  - can  be  accurately  accounted  for  as  described  below.  In  fact,  if 
the  integrand  could  be  represented  exactly  by  the  product  of  a simple  singular 
function,  such  as  and  a polynomial  of  order  2N-1,  then  using  N 

integration  points  would  produce  the  integral  with  no  error.  Although  the  true 
integrand  cannot  be  represented  exactly  in  the  above  manner,  by  taking  enough 
integration  points  the  error  may  be  mede  as  sm-.ll  as  desired. 

First,  the  pressure  influence  coefficient  of  Equation  (8.32)  is  rewritten 
by  performing  a change  of  variables: 


Cj  (3c, y) 


A(_L 

Rn  ddx 


r m rl 

J 

-m  -t  - 


“ +£  /•*  _ 

-ikAx  Tc  dAx 
e . cos  -JJ-  • ady. 


m ft  +2  _ 


r* 

ik 

•un  J 


e‘ik  Ax  . cos  ^ 

rf  R 


4 -ik  (i+i)/“+i 

‘ ^ 6 I 

■'m 


k R 

IT 


\ . MZ 


+ _i  U -i) 
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fm  +4 
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(8.34) 
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! 1 i I Receiving 

(a)  Type  of  Singularity  Point 

and  Integration  Area 1 j 

a None  - Unit  Square  I I 

b Along  2 Sides  - 1/1*  Square  d At  2 Corners  - Unit  Square 

c Along  1 Side  - 1/2  Square  e At  1 Corner  - Unit  Square 


where 


Ax  = x - § 

Ay  * y - i) 


£,m  are  the  components  of  the  distance  between  the  j 
point  (x,  y)  - see  Figure  8.4b. 

*^(i  - i)2  - (7  - ^i)2 


.th 


box  center  and  the 


Hi  =*  Jlii  - if  - (y  - ^)2 


4 

Next,  the  two  single  and  one  double  integrations  are  performed  using  various 
quadrature  formulae  of  Reference  8-9.  Referring  to  Figure  8. ha,  five  cases 
can  arise  depending  on  how  the  sending  box  is  cut  by  the  inverse  Mach  cone 
from  (x,  y): 

(a)  Box  not  touched:  L a 2,  m s l - 2. 

(b)  Box  at  apex  of  the  cone:  l 3 0,  m = 0. 

(c)  Box  split  by  the  cone:  / 3 m,  m > 0. 

(d)  Box  touched  at  two  corners  by  the  cone:  l 3 1,  m 3 0. 

(e)  Box  touched  at  one  comer  by  the  cone:  t > 1,  m 3 i - 1. 

In  case  (a),  since  there  is  no  singularity,  the  following  quadrature 
formula  (Equation  25.4.30,  Reference  8-9)  may  be  used  for  both  the  spanwise 
and  chordwise  integrations: 


J-b  n 

f(y)  dy  3 ^H.  « f(yt), 

a i=l 


where 


b - a 

yi  = T Xi 


b + a 


i,t 

xi  is  the  i zero  of  the  Legendre  polynomial,  Pn 
H , the  weighting  function  at  the  i quadrature  point, 

■ 2/d  - *,)  [p;  (x,)] 2. 


(8.35) 
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where 


yi  * a ♦ (b-a)  (l-x2) 

x^  is  the  i pos.  zero  of  the  Legendre  polynomial,  P2n< 

(2n) 

' are  the  Gaussian  weights  of  order  2n. 

Again,  Formula  (8.35)  is  used  in  the  chordwise  integration. 

For  case  (d),  an  integration  is  first  pevfoimed  over  the  triangular 
shaded  area  in  Figure  8.4a  (comprising  boxes  b,  c,  and  d)  using  the  quadra- 
ture formulae  for  case  (b).  Then  by  subtracting  the  previously  derived 
pressure  influence  coefficients  for  boxes  b and  c,  the  desired  integral  over 
d is  achieved. 

Finally,  for  case(e),  the  integral  for  an  aggregate  area  (see  dotted  area 
in  Figure  8.4a)  consisting  of  a triangle  on  the  Mach  line  and  the  subject  area 
is  determined;  and  from  it  the  integral  for  a triangle  is  subtracted.  When  the 
spanwise  integrations  are  performed,  Equation  (8.37)  is  employed,  while 
Equation  (8-35)  is  used  in  the  chordwise  direction. 

For  case  (a)  and  all  chordwise  integrations,  the  present  program  uses  six 
integration  points  for  the  quadrature.  In  the  spanwise  integrations,  six  points 
are  used  when  k < ^ (®j)2,  while  twelve  points  are  employed  when  k 2 j (| )2. 

At  a given  Mach  number  and  reduced  frequency,  the  pressure  influence 
coefficients  are  functions  of  only  l and  m - the  separation  between  the  sending 
and  receiving  box  centers.  Consequently,  influence  coefficients  are  computed 
by  the  above  formulae  only  once  for  each  admissable  i,  m pair  ( f 2 0, 
mil)  and  are  used  repeatedly  where  needed. 

The  pressure  on  any  box  is  a function  of  the  downwash  of  only  those  boyes 
within  the  inverse  Mach  cone  emanating  from  its  box  center.  For  a surface, 
the  edges  of  which  are  all  supersonic,  the  pressure  is,  furthermore,  only 
influenced  by  boxes  on  the  planform..  If  any  of  the  surface  edges  are  subsonic, 
however,  there  are  regions  adjacent  to  these  edges  which  do  affect  the  pressure 
of  some  areas  on  the  planform.  To  account  for  this  effect,  the  concept 
(Reference  8-10)  of  a permeable  "diaphragm"  is  introduced  in  these  regions. 

This  permeable  sheet  does  not  alter  the  flow  nor  can  it  sustain  pressure.  It 
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is  bounded  by  the  surface  edge  and  the  Mach  lines  emanating  from  the  corners 
of  the  lifting  surface  - see  Figure  8.5. 


The  relationship  between  the  pressure 
the  lifting  surface  and  the  diaphragm  area 


and  downwash  on  the  combination  of 
can  be  written: 


where 


pw 

fwW^CWD 

CDW|  CDD 

i 

1 • 


(8.38) 


Pw  is  the  pressure  on  the  wing  boxes 

PD  is  the  pressure  on  diaphragm  boxes 

are  the  influence  coefficients  giving  pressures  on  the  wing  boxes 
due  to  downwash  on  the  wing  boxes 

Cyp  are  the  influence  coefficients  giving  pressures  on  the  wing  boxes 
due  to  downwash  on  the  diaphragm  boxes 


CDW  are  the  influence  coefficients  giving  pressures  or.  the  diaphragm 
boxes  due  to  downwash  on  the  wing  boxes 


Cpjj  are  influence  coefficients  giving  pressures  on  the  diaphragm  boxes 
due  to  downwash  on  the  diaphragm  boxes 

ws  is  the  known  downwash  on  the  wing  boxes 

wD  is  the  unknown  downwash  on  the  diaphragm  boxes. 


Since  the  pressure  on  any  diaphragm  box  is  zero,  then 


[CBw]{Ws} + [CDd] {Wd|  * |0}'  (3'39) 

and  the  unknown  diaphragm  downwash  can  be  evaluated  by 

h(-  * [cdd]  UM  <8-Uo> 
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Substituting  this  result  into  Equation  (8.38),  yields  the  final  expression  for 
the  pressures  on  the  wing  boxes: 


where 


M ■ MM 

H ■ M - M M’1  [cd»] 


(8.41) 

(8.42) 


For  maximum  computer  efficiency,  the  actual  calculation  of  the  pressures, 
Pw,  is  generally  performed  in  a different  manner  from  that  implied  by  Equation 
(8.4i).  The  calculation  of  the  aerodynamic  influence  coefficient  matrix,  ^Aicj 
leads  to  either  extensive  use  of  core  in  storing  matrices  or  a large  number  of 
I/O  operations  if  the  matrices  are  stored  on  data  devices.  If  there  is  no  neeu 
for  saving  the  [mcJ  array  for  future  use,  the  machine  operations  can  be 
appreciably  reduced  by  computing  the  pressures  as  follows: 
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In  this  way,  core  storage  requirements  are  minimized,  since  only  a vector  need 
be  stored  in  going  from  one  step  to  another. 
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Cnee  the  oscillatory  pressures  are  computed,  the  normalized  generalized 
aerodynamic  forces  are  computed  by 


hr(x,  y)  C (x,y)  dx  dy. 


(8.U4) 


Since  the  boxes  are  assumed  to  be  very  small,  the  integration  can  be  replaced 
by  a summation  over  every  box  on  the  planform: 
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(8.45) 


vhere 

h^  is  the  deflection  of  the  j*'1’  box  in  the  rth  mode 

is  the  pressure  coefficient  of  the  j^*1  box  in  the  sth  mode 
is  the  area  of  the  jth  surface  box. 

It  should  be  noted  that  each  box  area  is  b /g.  The  program  automatically 
establishes  the  gridwork  of  boxes:  From  a user  specified  number  of  boxes 

desired,  the  program  calculates  the  box  size  necessary  for  the  boxes  to  cover 
the  planform  and  diaphragm  and  to  align  with  the  inboard  and  outboard  planform 
edges.  Consequently,  no  boxes  overhang  the  planform  side-edges.  (Referral 
to  Figure  8.5  here  and  in  the  remaining  discussion  may  be  helpful.)  Each  box 
is  designated  a3  either  a wing  box  or  diaphragm  box  depending  on  whether  its 
center  lies  on  or  off  the  wing,  respectively.  In  general,  boxes  dfi  overhang 
the  leading  and  trailing  planform  edges  causing,  in  effect,  a jagged  represent- 
ation of  these  edges.  For  most  configurations,  this  jaggedness  has  been  found 
to  have  little  effect  on  the  accuracy  of  the  computed  generalized  aerodynamic 
forces,  providing  that  the  box  grid  is  not  too  coarse.  It  follows  that 


since  pressure  is  assumed  constant  over  each  box,  wing  boxes  that  overlap 
the  planfora  edges  support  a force  on  their  off-wing  portion  as  well  as  the 
portion  actually  lying  on  the  surface.  This  results  in  a computed  force  in 
excess  of  the  correct  amount.  However,  diaphragm  boxes  that  overlap  the 
planform  edges  support  no  force  on  either  the  on-  or  off-wing  portions  and, 
hence,  result  in  a deficit  in  the  computed  wing  forces  that  tends  to  balance 
the  excess  obtained  from  overhanging  wing  boxes.  In  Figure  8.5,  this  balance 
is  illustrated  for  one  representative  pair  of  boxes.  The  shaded  area  is  the 
off -wing  portion  of  wir.g  box  A while  the  crosshatched  area  is  the  roughly 
balancing  on-wing  portion  of  diaphragm  box  B. 

When  a highly  swept  surface  is  analyzed  for  a relatively  low  Mach  number, 
the  number  of  forward  diaphragm  boxes  can  become  so  great  as  to  cause  an 
appreciable  increase  in  computing  time.  However,  the  downwash  in  the  diaphragm 
region  decreases  very  rapidly  in  the  forward  streanwise  direction.  To  save 
computing  time,  the  present  program  takes  this  rapid  decay  of  diaphragm  down- 
wash  into  account  and  allows  the  user  to  request  a box-elimination  option 
whereby  the  diaphragm  boxes  are  ignored  forward  of  a user  specified  distance 
ahead  of  the  leading  edge. 

In  the  program,  provision  is  made  for  computing  aerodynamic  force  co- 
efficients and  center-of-pressure  locations.  The  user  may  use  this  facility 
to  compare  known  steady  state  data  with  computed  values  to  determine  the  number 
of  boxes  required  for  a satisfactory  sol  tion.  Another  approach  is  to  vary  the 
number  of  boxes  and  look  for  convergence  in  the  stability  coefficients. 

8.2.3  Subsonic  Doublet-lattice  Program 

For  the  added  capability  in  the  subsonic  regime  of  analyzing  control  sur- 
face configurations,  multiple  interfering  surfaces  and  interfering  surface-body 
configurations,  the  doublet -lattice  program  of  Reference  8-11  is  used.  The  formula- 
tion of  this  method  is  different  from  the  assumed-pressure-function  method  but 
starts  with  sane  integral  equation  relating  the  wash  normal  to  a harmonically 
oscillating  surface  to  the  lifting  pressure: 

v (Pj)  * ^ JJ  Cp<P)  K^,  P,  k,  M)  d?  do  , (8.2) 
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M*--*  T.' 


where 


P is  any  point  on  the  planform,  the  coordinates  of  which  are  §,  T|,  £ 

P is  a j point,  the  coordinates  of  which  are  x,  y,  z 

x,  5 are  streamwise  coordinates 

y,  1)  are  spanwise  coordinates 

z,  £ are  vertical  coordinates 


a 

w 


is  the  tangential  spanwise  coordinate 


w It 

is  the  normalwash  angle  *u  = rt  + *b  ^ 

o 


is  the  differential  pressure  coefficient 


is  the  kernel  function  relating  the  normalwash  at  Pj  to  the  a unit 


h 

a 

k 

S 


pressure  at  P 

is  the  Mach  number 

is  the  reference  semichord 

is  the  vertical  displacement  of  the  surface 

is  the  streamwise  slope  of  the  deformed  surface 

is  the  reduced  frequency 

is  the  surface  area  of  all  lifting  surfaces  included  in  the  analysis. 


If  the  surface  is  divided  into  J elements  over  which  the  pressure  is  assumed 
constant,  the  previous  equation  becomes 


J 

w (x,  y,  z)  = JJk  (x-S,  y-D,  z-C,  to,  M)  d§  do. 

3 1 T.  I W.IW.1W1 


(8.46) 
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NexJ  , the  pressure  is  assumed  to  arise  from  a loaded  line  at  the  i-chord 
of  each  element.  As  illustrated  in  Figure  8.6,  this  is  equivalent  to  an  un- 
steady horseshoe  vortex  whose  hound  portion  lies  along  the  chord  of  the 
element.  The  resulting  expression  is: 

J 

C / K <x"Sl*  y-1'*  Z-C>  «,  M)  do,  (8.47) 
pj  3 J 

ELEMENT 
i 

where  4^  is  the  length  of  the  average  chord  of  element  j,  and  the  integration 
is  taken  along  the  £ -chord  line  of  the  element. 


v (x,  y,  z)  - 

J-l 


Figure  8.6.  Horseshoe  Vortex  Element  Used  in 
Doublet- Lattice  Method. 
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If  ':he  downwash  is  calculated  at  points  located  at  the  3A-chord  midspan 
of  each  of  the  J elements  and  the  above  equation  is  satisfied  at  each  of  these 
points,  the  follc-'ing  set  of  equations  holds  for  i * 1 to  J: 

J 


i 

ELEMENT 

d 


This  can  be  written  in  matrix  form  as: 

H - m {cp}>  (8.49) 

where  relates  the  dcwnwash  at  the  A*1  point  to  the  pressure  over  the 
element.  By  solving  this  set  of  linear  equations,  the  pressure  distribution 
over  the  surface  is  calculated. 

Because  the  doublet  lattice  method  is  not  explicitly  fitting  pressure 
functions  on  a planform  as  is  done  in  the  assumed-pressure-function  method, 
multiple  aerodynamically  interacting  surfaces  can  be  modelled  by  simply  de- 
fining lattice  elements  on  each  surface.  In  this  case,  Equation  (8.48)  still 
holds,  but  J is  now  the  number  of  elements  on  all  surfaces. 

To  account  for  aerodynamic  interaction  between  bodies  and  surfaces, 
Reference  8-11  describes  a method  of  modelling  each  body  by  axial  doublets 
along  the  body  axis  and  by  panels  of  unsteady  horseshoe  vortex  elements  on 
the  body  surface  in  the  vicinity  of  each  lifting  surface  with  which  the  body 
might  interact  (see  Figure  8.7,  for  example).  The  strength  of  each  axial 
doublet  is  calculated  by  slender  body  theory.  The  incremental  downwash  on 
the  panels  on  the  body  surfaces  and  on  the  lifting  surfaces  caused  by  these 
axial  doublets  is  then  computed  and  subtracted  from  the  prescribed  downwash 
for  these  surfaces.  Equation  (8.49)  becomes 
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'■TO1 


where 


p(S)  is  the  pressure  coefficient  distribution  on  the  liftiig  surfaces 


P(I) 


P(B) 


LSB 


IB 


SS 


SI 


"IS 


"II 


Since 


is  the  pressure  coefficient  distribution  on  the  interacting  body 
surfaces 

is  the  pressure  coefficient  distribution  on  the  slend-’r-body 
axial  elements 

is  the  downwash  distribution  prescribed  on  the  lifting  surfaces 

is  the  downwash  distribution  prescribed  on  the  interacting  tody 
surfaces 

i3  the  downwash  distribution  on  the  lifting  surfaces  caused  by  unit 
pressure  coefficients  along  the  slender  bodies 

is  the  downwash  distribution  on  the  interacting  body  surfaces  caused 
by  unit  pressure  coefficients  along  the  slender  bodies 

r 

is  the  downwash  distribution  on  the  lifting  surfaces  caused  by  unit 
pressure  coefficients  on  the  lifting  surfaces 

is  the  downwash  distribution  on  the  lifting  surfaces  caused  by  unit 
pressure  coefficients  on  the  interacting  body  surfaces 

is  the  downwash  distribution  on  the  interacting  body  surfaces  caused 
by  unit  pressure  coefficients  on  the  lifting  surfaces 

is  the  downwash  distribution  on  the  interacting  body  surfaces  caused 
by  unit  pressure  coefficient  on  the  interacting  body  surfaces. 


P(B) 


was  calculated  from  the  downwash  and  geometry  of  each  body  using 


slender-body  theory,  the  only  unknowns  in  this  matrix  equation  are  { C 


P(I) 


P(S) 


and 


This  set  of  equations  is  solved  for  these  pressure  distributions  by 


a standard  linear  system  solution  algorithm  using  Gaussian  triangularization 
and  back  solution  (see  Reference  8-3,  pp  428-429). 


As  in  the  Mach-box  method,  generalized  aerodynamic  forces  arising  from 
lifting  surface  pressures  are  calculated  by 


LhrJ'  C-"  A' 


rs 


Si)  P(J)  *(J) 


(8.U5) 


J-l 


where,  here,  is  the  area  of  the  element.  This  can  be  rewritten  in 

matrix  form  for  all  N modes  as: 


(8.51) 


(NXN) 


(NXJ)  (JXJ)  (JXN) 


where  JSyJ  Is  a diagonal  matrix.  When  interacting  bodies  are  present,  three 
generalized  forces  are  present: 


[«]  ‘ [«5  * 5I  * 


(8.52) 


The  first  arises  from  the  pmiuct  of  deformations,  pressures  and  areas  of  the 
lifting  surface  elements;  the  second  from  the  product  of  those  of  the  inter- 
acting body  surface  elements;  and  the  third  from  the  product  of  those  of  the 
body  elements  - the  appropriate  areas  in  this  last  case  are  the  products  of  the 
body  element  diameters  and  le.tgths. 


To  obtain  satisfactory  pressure  distributions,  the  lifting  surface  must 
be  divided  into  strips  of  elements  whose  edges  are  parallel  to  the  free  stream. 
An  example  is  shown  in  Figure  8.8.  Additionally,  element  edges  should  lie  along 
surface  edges,  fold  lines  and  control  surface  hinge  lines.  Three  guidelines 
should  be  observed  in  subdivision: 


(1)  The  leading  and  trailing  edges  of  adjacent  pairs  of  elements  should 
be  aligned  and  located  at  a constant  percent  of  the  strip  chord  when 
possible. 

(2)  The  dimensions  of  elements  should  be  decreased  in  the  directions  and 
regions  of  large  gradients  in  pressure  and/or  downwash,  such  as  near 
hinge  lines,  leading  edges  and  wing  tips. 
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Figure  8.8.  Doublet-Lattice  Modeling  of  a Delta  Wing  with 
Cranked'  Leading  Edge  and  Two  Control  Surfaces. 

(3)  The  aspect  ratio  of  each  element  should  be  unity  or  less.  However, 
this  is  not  always  possible, especially  in  regions  where  a large 
pressure  gradient  is  expected.  This  is  evident  in  the  example 
shown  in  Figure  8.7. 

The  optimum  configuration  is  predicated  by  the  need  for  keeping  the  number 
of  elements  to  a minimum  and  at  the  same  time  generating  generalized  air  force 
terms  that  are  satisfactory  for  flutter  analysis.  It  may  be  necessary  to  test 
a number  of  trial  configurations  and  compare  results  before  making  a final 
analysis  with  a fixed  element  layout. 

8.2.4  Modal  Interpolation 

For  each  of  the  above  aerodynamics  routines,  the  normalwash  angle  is 
required  at  specified  aerodynamic  grid  points  on  the  lilting  surface: 
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ij 


Wij 

U ~ "ij 


(8.53) 
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where  a.  . arid  h, , are  the  streamwise  slope  and  the  displacement  of  point  i in  I 

the  y vibration  mode.  Since  the  modal  data  from  the  vibration  analysis  is 
specified  on  a dynamics  grid  which,  in  general,  does  not  coincide  with  the  aero-  J 

dynamics  grid,  an  interpolation  is  required.  The  procedure  used  consists  of 
representing  the  dynamics  grid  as  a set  of  sp&nwise  oriented  lines  connecting 
nnde  points  at  which  modal  deflections  are  specified.  These  deflections  are 
interpolated  along  the  lines  to  each  spar.wise  station  at  which  aerodynamic 
grid  points  lie  and  then  are  interpolated  chordwise  along  each  such  station  to 
each  aerodynamic  grid  point.  This  scheme  is  illustrated  in  Figure  3.9. 


Figure  3.9*  Modal  Interpolation  Scheme. 

In  a variation  of  this  scheme,  available  as  a pregram  option,  modal 
streamwise  slopes,  as  well  as  deflections,  are  specified  as  input  data  along 
a single  spanwise  line.  The  program  then  creates  a second  line  parallel  to 
and  at  a specified  streamwise  distance  from  the  specified  line  and  transforms 
the  modal  slopes  and  deflections  to  a set  of  deflections  alone  each  line: 


where 


| are  tne  deflections  and  stresmwise  slopes  prescribed 
along  the  first  line 

are  the  deflections  computed  along  the  second  line 
is  an  identity  matrix 

is  a diagonal  matrix  having  the  specified  separation  between 
the  lines  as  each  of  its  diagonal  terms. 


Using  this  set  of  lines  and  deflections,  the  interpolation  proceeds  according 
to  the  original  scheme. 

for  a surface  with  a control  surface  or  some  other  special  region  across 
the  boundaries  .of,  which  the  modal  deflections  are  r.ot  continuous,  an  option  i3 
available  whereby  the  interpolation  is  performed  over  the  main  surface  and  over 
the  control  surface(s)  separately  u*.ing  separate  sets  of  lines  and  node  points. 

In  this  manner,  modal  discontinuities  across  the  boundaries  are  preserved. 

The  calculations  performed  in  the  above  schemes  use  the  Lagrangian  inter- 
polation formula  of  Reference  8-3,  pages  60-68.  Accordingly,  a polynomial,  g(x), 
is  determined  as  an  approximation  to  a function,  f(x),  the  value  of  which  is 
known  at  each  of  Ii+1  points,  |xq,  polynomial  is  computed 

by 


g(x) 


(x“xo} 


k*0 


(x.-x  j 
* o 


1 Oc-Vl)  U-Vp)  •••  (x-x:;) 
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f(xk). 


(8.55) 


In  the  present  program,  the  polynomial,  g(x),  is  limited  to  degree  '.I  = 3 to 
minimize  convolutions  in  the  approximation.  Consequently,  ir.  cases  for  which 
the  modal  deformations  are  specified  at  more  than  four  points  spanwise  or  chord- 
wise,  the  approximation  is  a pier evise-continuous  cubic  polynomial.  ■•Is  a program 
option,  the  user  can  further  restrict  the  polyr.cxr.ial  to  a linear  or  parabolic 
function  in  regions  where  extrapolation  is  needed  either  spanwise  or  chordwise. 


8.2.5  .Multiple  Surface  Capabilities 

Ir.  the  doublet-lattice  procedure,  multiple  aerodyr.amicaily  interacting 
surfaces  may  be  treated  as  described  in  Section  8.2. 3.  Although  the  Mach  box 
procedure  and  the  assuned-pressure-function  procedures  do  not  account  for 
interaction  effects,  these  methods  may  be  used  to  analyse  a multiple  surface 
vehicle  where  aerodynamic  interaction  ic  not  present.  This  is  achieved  by 
calculating  the  generalised  aerodynamic  force  contributions  for  each  of  the 
vehicle  lifting  surfaces  and  adding  these  together  before  passing  them  to  the 
flutter  solution  routine.  Additionally,  the  user  can  request  that  the  aero- 
dynamic forces  of  a particular  surface  be  multiplied  by  a specified  scalar 
before  the  addition  is  made,  so  that  empirical  knowledge  of  the  actual  lift 
of  a surface  as  compared  to  its  calculated  value  may  be  reflected. 

8.3  GENERALIZED  AERODYNAMIC  FORCE  INTERPOLATION 

For  a flutter  analysis,  generalised  aerodynamic  forces  are  required  at 
several  reduced  velocities.  To  reduce  the  computation  time  needed  to  obtain 
these  forces,  interpolation  is  used  to  determine  the  generalised  aerodynamic 
forces  at  all  desired  reduced  velocities  from  forces  directly  computed  at  a 
small  set  of  selected,  reference  reduced  velocities.  A separate  interpolation 
is  performed  on  the  real  and  on  the  imaginary  part  of  each  term  of  the  general- 
ised aerodynamic  force  matrix.  As  in  the  modal  interpolation  of  section  8.2.L, 
lAgrangian  interpolation  {Reference  8-3>  pages  60-68)  is  used  and  the  approxi- 
mating function  is  limited  to  a piecewise  continuous  second-  or  third-order 
polynomial, 

In  the  present  program,  the  user  supplies  a goodness-of-fit  tolerance 
and  six  reference  reduced  velocities  ordered  by  increasing  value  and  distribute 
over  the  range  required  in  the  subsequent  flutter  analysis:  Jv^,  v^,  ...  Vg  J 

dsing  generalised  forces  computed  at  three  of  these  six,  a generalized  aero- 
dynamic matrix  at  a fourth  is  determined  by  interpolation  and  compared  with  a 
matrix  of  generalized  forces  directly  computed  at  that  reduced  velocity.  The 
comparisons  are  performed  for  various  arrays,  C,  formed  from  combinations  of 
the  terms  of  the  generalized  force  matrix  and  take  the  form  of  the  test: 


After  the  goodnesa-of-fit  tests  have  been  made,  the  procedure  for 
obtaining  generalized  forces  at  any  selected  reduced  velocity  is  summarized 
in  Table  8.3. 


TABLE  5.3.  CEERALIZED  AERODYNAMIC  FORCE  INTERPOLATION 
TROCEDURE  IN  VARIOUS  REDUCED  VELOCITY  RANGES 


TEST  A 
PASSED 


TEST  B 
PASSED 


TEST  C 

TEST  C 

PASSED 

FAILED 

Extrapolation 
using  parabolic 
fit  through 

^1»  ®3» 


Interpolation  Interpolation  Interpolation 

using  parabolic  using  cubic  using  cubic 

fit  through  fit  through  fit  through 

^3>  ^6  S2»  ^3>  ^6  §l>  ^3,  Q4 


Interpolation 
using  cubic 
fit  through  _ 

\>  ^3*  V \ 


Interpolation 
using  cubic 
fit  through  _ 

\>  \>  V \ 


Interpolation 
using  cubic 
fit  through  _ 

^3»  %>  % 


fit  through 
^3>  ^6 


Extrapolation 

Extrapolation 

Extrapolation 

using  parabolic 

using  parabolic 

using  parabolic 

fit  through 

fit  through 

fit  through 

\t  §3>  Qg 

^3»  ^4*  ^6 

^5>  *<6 

Note:  v^,  * (l/k)^,  ordered  by  increasing  reduced  velocity 


If.  choosing  the  reference  reduced  velocities,  the  following  considerations 
should  be  observed: 

(1)  For  the  computed  aerodynamic  forces  to  be  most  accurate,  the  lowest 
reduced  velocity  should  as  high  as  possible  consistent  with  the 
range  of  values  needed  for  the  flutter  solutions. 

(2)  Generally,  the  first  and  sixth  reference  reduced  velocities  should 
span  the  expected  range  of  v for  which  generalized  aero  forces  are  to 
be  calculated.  This  precaution  eliminates  the  need  for  extrapolation. 

(3)  The  reference  reduced  velocities  should  emphasize  regions  where 
flutter  is  expected. 

8.4  FLUTTER  SOLUTION  PROCEDURES 

There  are  two  solution  options  available  to  the  user  for  the  flutter 
analysis: 

(1)  the  conventional  k-method,  using  the  QR  algorithm  (Reference  8-12, 
pages  515-568)  to  determine  eigenvalues. 

(2)  the  p-k-method,  using  a determinant  iteration  procedure  with  a 
quadratic  predictor  (Reference  8-12,  page  435)  to  determine 
eigenvalues. 

Generalized  aerodynamic  forces  required  at  various  reduced  velocities  in  these 
analyses  are  usually  determined  by  the  interpolation  procedure  described  in 
Section  8.3.  However,  if  the  user  desires,  directly  computed  aerodynamic 
forces  can  be  used  in  the  k-analysis.  To  help  the  analyst  study  the  flutter 
mechanism,  various  options  such  as  the  automatic  preparation  of  plots  of  modal 
structural  damping  and  frequency  as  functions  of  airspeed  are  available. 

If  the  redesign  capability  of  FASTOP  is  to  be  used,  the  p-k  flutter 
analysis  must  be  selected,  in  which  case  various  items  such  as  the  flutter 
vectors  are  computed. 
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8.4.1  Flutter  Equation 


The  basic  differential  equation  for  a lifting  surface  in  oscillatory 
motion  is: 

[m]  [q  j " [<*]  |<l  | + [K]  {<1  ) - O.  (8.57) 

where 

is  the  generalized  mass  matrix, 
is  the  generalized  aerodynamic  force  matrix, 
is  the  generalized  stiffness  matrix, 
is  the  generalized  coordinate  vector. 

The  generalized  aerodynamic  force  matrix  equals  the  product  of  the  aero- 
dynamic matrix,  £aJ,  and  the  dynamic  pressure,  ipU2,  where  an  element  of  £aJ 
is  defined  as: 
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(8.58) 


where  hr(*x,  y)  is  tre  deflection  diatritxition  in  the  rth  mode,  C (x,  y) 

th  P(s) 

is  the  pressure  coefficient  distribution  in  the  s mode,  and  the  integration 
is  performed  over  the  lifting  surface  area,  S.  Hence,  Equation  (8.57)  can  be 
rewritten  as 


H H -tH  M *M  M*°-  (s-59) 

If  structural  damping  is  present  in  the  system,  the  generalized  stiffness 
array  should  be  modified  as: 

[k]  - [l  + ic]  [k],  (8.60) 
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where  £gJ  is  a diagonal  matrix  of  modal  damping  and  £lj  is  the  identity  matrix. 


The  stiffness  in  the  j mode  is  then 


’ C1  * 1 «.(!))  V 


(a. 61) 


where  is  the  hysteretic  or  structural  damping  in  the  j " mode.  Then 

Equation  (8.59)  becomes 


M M-M  M-o. 


(8.62) 


For  harmonic  motion,  the  generalized  coordinates  may  be  written  as 


q = ( q I eAt,  where  q is  time  independent  and  a may  contain  damping. 


Equation  (8.62)  can  then  be  rewritten  as: 


[>  W -#[»]*  w]  I 


q I * o. 


(8.63) 


8.  **.2  k-Method 


In  the  "American  approach"  (or  k-method),  since  the  conventional  aero- 
dynamic theories  are  valid  only  for  undamped  oscillations,  the  aerodynamic 
matrix  is  computed  for  a chosen  reduced  frequency,  k,  and  the  a is  considered 
to  be  undamped  (a  * iou): 


p [m]  - ioU2  [A(k>]  ♦ [k]  j 


q * 0. 


(8.64) 


This  equation  may  be  rewritten  as 
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where  the  reduced  frequency  is  defined  as 


k-SS,  (8.66) 

where 

b is  a referer*'  senichord 

to  is  the  freer.  ;ncy  of  oscillation 

U is  the  velocity  of  the  free  stream. 


Now.  an  unknown  hysteretic  structural  damping  term,  g,  is  introduced  to 
bring  the  system  into  neutral  stability.  Applying  this  damping  to  the  jth 
diagonal  term  of  the  stiffness  matrix  of  Equation  (8.65)  results  in: 


* 2 (-1  + igs(j)  + lg) 

V 

(8.67) 

For  small  values  of  damping  ^g  • gs(j)<<3)>  this 
be  approximated  by 

on-diagonal  term  may 

KJj* 

(8.68) 

* x (1  * igs(d))  Kjj  “ x 

where 

•V 

(8.69) 

X * -|r  (1  + ig). 

(8.70) 

ui 


Substitution  in  Equation  (8.65)  yields 

[l  M - M • H?  M] 


Defining  the  aerodynamic  term  as.  Q (k),  this  becanes 


Solving  this  equation  for  each  eigenvalue,  one  obtains  the  values  of  the 
unknown  damping  and  the  accompanying  u)  of  the  system.  A positive  vel  je  of 
the  damping  means  that  structural  damping  must  be  added  to  the  system  to  bi ing 
it  to  neutral  stability;  i.e. , the  system  is  unstable.  A negative  value  has 
the  opposite  interpretation. 

A flutter  analysis  by  the  k-raethod  begins  with  the  calculation  (or  inter- 
polation) of  generalized  aerodynamic  forces,  [5(k)]  , for  a chosen  set  of  reduced 
frequencies.  Using  the  QR  algorithn  of  Reference  8-12,  pages  515-568,  eigen- 
values are  determined  from  Equation  (8.72)  for  each  reduced  frequency.  Using 
Equations  (8.70)  and  (8.66),  the  reduced  frequency  and  the  real  and  imaginary 
parts  of  each  eigenvalue  are  converted  to  a frequency,  a required  damping,  and 
a free  stream  velocity.  9y  plotting  these  dampings  and  frequencies  as  functions 
of  these  velocities,  the  critical  flutter  speed  (the  lowest  nonzero  velocity  at 
which  g » 0)  and  the  accompanying  flutter  frequency  are  determined. 

Having  derived  the  eigenvalues  of  Equation  (8.72),  the  eigenvectors  are 
determined  using  ar.  inverse  iteration  procedure  (Reference  8-12,  pages  619-625). 
These  are  used  to  establish  the  uniqueness  of  the  eigenvalues  in  a routine  using 
the  Gershgorir.  theorem  (Reference  8-12,  pages  638-6U6). 

3ecause  of  the  assumptions  implicit  in  this  approach,  the  subcritic* 1 damp- 
ing and  frequency  trends  are  generally  inaccurate.  Occasionally,  the  methou 
produces  a multiple  valued  function  of  damping  vs.  velocity,  making  it  difficult 
to  order  the  roots  in  a routine  to  automatically  determine  the  flutter  speed. 

The  advantage  of  this  approach,  however,  is  its  speed;  solutions  to  linear 
eigenvalue  problems  are  relatively  easy  to  compute. 

S.h.3  p-k  Method 

An  alternate  approach  to  the  solution  of  the  flutter  equation,  which  gives 
better  subcritical  trends  and  does  not  lead  to  double  valued  functions  of  damp- 
ing vs.  velocity,  is  the  p-k  method  of  Reference  8-1.  The  generalized  coor- 
dinates are  assumed  to  be  damped  harmonic  functions;  hence,  a becomes 


a * a *•  iuj 
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where  v is  the  darling  coefficient  equal  to  — in  — — (with  a and  a , as  the 

c<«  a n r.-i 

n 

amplitudes  of  successive  cycles  of  oscillation).  Equation  (8.63)  then  becomes 

where  the  available  constant  amplitude  (undamped)  aerodynamic  theories  must 
be  used  to  compute  the  aerodynamic  matrix. 

Unlike  the  k-nethod,  an  air  speed  U is  row  selected  for  which  Equation 
(S.7M  will  be  solved.  The  characteristic  equation  for  the  eigenvalues , p, 
is  then  written  as 


F<P.  «*)  * | (j*.f  P2  [m]  * [k]  - [a( k )]  j = 0.  (9.75) 

In  Reference  8-1,  this  equation  is  solved  for  p by  the  iterative  Regula  Falsi 
method.  However,  it  is  noted  in  that  reference  that  this  algor  it  hr.,  which 
uses  a linear  predictor,  occasionally  exhibits  nonconverger.ee.  Consequently, 
in  the  present  program  this  algorithm  is  replaced  by  Muller 1 s method 
(Reference  6-12,  pages  •+35-U3S),  which  employs  a quadratic  predictor.  This 
means  that  each  root  estimate  is  a function  of  the  most  recent  three  estimates. 

In  preparing  data  for  the  p-k  solution  option,  the  user  - by  specifying 
ar.  initial  eirspeed,  an  airspeed  increment  ar.d  the  total  number  of  airspeeds 
for  which  solutions  are  desired  - defines  the  velocity  range  over  which  the 
flutter  solution  is  to  be  determined.  At  the  initial  velocity,  the  program 
makes  three  estimates  of  each  root,  calculates  corresponding  generalized 
aerodynamic  forces  ar.d  computes  F,  the  flutter  determir.ar.t  of  Equation  '5.75), 
for  each  estimate.  After  fitting  a quadratic  to  the  three  F's,  the  program 
determines  the  value  ox’  the  root  for  which  this  quadratic  equals  ;ero.  This 
value  becomes  the  next  estimate  of  the  root.  Using  the  three  most  recent  root 
estimates,  this  iteration  is  continued,  determining  new  estimates,  ■.■.‘her.  the 
iteration,  has  converged,  the  process  is  concluded  for  this  root  and  begun  for 
the  next. 
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Solutions  are  determined  for  each  specified  velocity  in  the  user-defined 
range,  using  roots  determined  at  the  previous  velocity  as  initial  estimates 
at  the  next  velocity.  Optionally,  the  user  can  instruct  the  program  to  obtain 
the  initial  estimates  by  extrapolation  from  roots  at  the  previous  two  velocities. 
Whenever  significant  changes  in  the  slope  of  the  damping  or  frequency  are 
detected  over  the  usei  -defined  interval,  the  program  automatically  obtains 
solutions  at  finer  velocity  increments.  (See  Figure  8.10.)  Finally,  when  a 
sign  reversal  in  the  damping  part  of  the  root  occurs  between  consecutive 
velocities,  the  program  initiates  a search  to  determine  the  velocity  for  which 
the  damping  is  zero,  i.e.,  the  flutter  speed. 

8.4.4  Flutter  Analysis  Features  and  Options 

To  help  the  user  study  the  flutter  mechanism,  the  program  allows  for  the 
variation,  in  a single  computer  submission,  of  the  stiffness  of  a chosen  mode 
and  of  the  number  of  normal  modes  (from  an  original  set)  included  in  the 
analysis.  Furthermore,  the  flutter  results  - root  damping  and  frequency  as 
functions  of  airspeed  - are  presented  in  graphical  (print-plot)  form  as  an 
integral  part  of  the  computer  listing.  In  addition,  there  is  provision  for 
obtaining  CALCOMP  plots  of  the  flutter  solution. 

Another  option  is  the  provision  for  listing  modal  components  of  the  eigen- 
vectors (flutter  vectors)  for  a given  velocity  or  reduced  velocity  range 
(depending  upon  whether  the  p-k  or  k method  is  used  for  the  solution).  If 
requested  by  the  user,  the  equivalent  physical-coordinate  vectors  at  the 
locations  specified  by  input  vibration  data  are  also  calculated  and  displayed. 

Optionally,  the  user  can  make  changes  to  the  terms  of  the  flutter 
equation  by  adding  structural  damping  and  by  revising  the  generalized  inertial 
or  stiffness  matrices. 

A divergence  analysis  can  also  be  performed,  using  the  aerodynamic  forces 
derived  from  a non-oscillatory  condition.  Ir.  which  case,  Equation  (8.59)  becomes: 


'4  [a  (k-o)]  M*H  M -a- 


(8.76) 
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Both 


[a]  and  [k] 


are  real  matrices  and  the  eigenvalue  to  be  determined  by  the 

w -■  ».  J o 

OR  algorithm  is  ir  in  this  case.  If  no  positive  eigenvalues  exist,  i.e.,  if 
£ 0,  the  surface  does  not  diverge;  otherwise  the  divergence  speed  is  the 
square  root  of  the  lowest  positive  eigenvalue. 


8.4.5  Flutter  Vector < for  Redesign 

When  redesign  is  to  be  done  in  FASTOP,  the  flutter  vector  (among  other 
quantities)  is  required.  The  flutter  vector  is  the  eigenvector  of  Equation 
(8.74)  for  the  critical  root  at  the  flutter  speed  - the  lowest  non-zero  speed 
at  which  the  damping  (the  real  part  of  the  root)  is  zero.  As  described  in 
Section  8.4.3,  the  program  automatically  searches  to  determine  the  flutter 
speed.  The  eigenvalue  p,  of  the  critical  root  at  this  speed  is  determined 
and  the  corresponding  eigenvector,  {a  { (the  flutter  vector),  is  found  by  the 
inverse  iteration  technique  (Reference  8-12,  pages  619-625). 


By  transposing  the  flutter  matrix  of  Equation  (8.75)  and  determining 
another  eigenvector  for  p,  the  associated  row  vector,  \ v P>  is  obtained. 
Subsequently,  the  following  parameter  is  formed  for  use  in  the  computation  of 
the  flutter  derivatives  (see  Section  9): 
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(8.77) 


where  Q wes  defined  by  Equations  (8.?l)  and  (8.72).  The  derivative 

ok 

is  obtained  by  differentiation  of  the  generalized  aerodynamic  force  inter- 
polation polynomial,  described  in  Section  8.3. 
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Section  9 


RESIZING  FOR  COMBINED  FLUTTER  AND 
STRENOTH  REQUIREMENTS 


9.1  SUMMARY 

Several  methods  for  sizing  the  finite  elements  of  a lifting-surface  struc- 
tural idealization  to  achieve  minimum  weight  design  under  combined  strength  and 
flutter-speed  requirements  were  developed  and  evaluated.  Two  basic  categories 
were  considered:  methods  based  on  a combination  of  energy  principles  and  opti- 
mality criteria,  and  procedures  employing  numerical-search  techniques.  Drawing 
upon  the  experience  gained  from  studies  of  both  of  these  basic  methods,  a re- 
sizing algorithm  was  developed  that  employs  a uniform-flutter-velocity-deri>*tive 
optimality  criterion  for  .flutter-critical  elements  and  the  fully-stressed-deslgn 
criterion  for  strength-critical  elements.  The  final  result  is  a practical,  auto- 
mated approach  for  dealing  with  large-scale  idealizations  having  both  structural 
and  mass-balance  design  variables. 

This  sectlor  provides  0.  summary  of  the  major  findings  from  the  evaluation 
of  the  candidate  flutter  resizing  methods  and  the  factors  that  led  to  the  selection 
of  the  final  algorithm,  which  is  discussed  in  detail.  A more  complete  description 
of  the  methods  examined,  their  underlying  theory  and  assumptions,  and  the  results 
they  produced  for  a representative  example  wing  structure  are  presented  in  Refer- 
ence 9-1. 

9.2  EVALUATION  OF  FLUTTER  RESI-  TIG  ALGORITHMS 

For  a structure  subjected  to  a single  flutter-speed  constraint,  and  no 
other  constraints  such  as  those  imposed  by  strength  and  mir.imun-manufacturing- 
gage  requirements,  it  can  be  shown  that  for  a minimum-weight  design,  the  deriva- 
tives of  the  critical  flutter  speed,  V^.,  with  respect  to  the  design  variable 
weights,  must  be  equal;  that  is 

— -=•  * constant  (9.1) 

Smi 

for  all  "i"  structural  finite  elements  and  mass-balance  weights. 

This  "optimality  criterion"  provides  a standard  of  local  optimality  under 
the  limited  condition  of  a single  design  constraint.  However,  realistic  structural 
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design*  involve  many  constraints  in  combination  (e.g.,  strength,  flutter  speed, 
minimum  manufacturing  gages)  and  their  respective,  separately-governing  opti- 
mality criteria,  oust  be  blended  into  a composite  criterion.  Moreover,  to  be  of 
real  value,  it  oust  be  possible  to  satisfy  such  a composite  criterion  by  a practi- 
cal and  efficient  resizing  procedure.  These  considerations  have  strongly  influ- 
enced the  development  of  the  combined  flutter/strength  resizing  algorithm  in 
FASTOP. 

The  finally  selected  approach  evolved  from  an  extensive  study  of  two 
classes  of  optimization  methods.  The  first  included  optimality-criteria  methods 
based  on  energy  concepts;  the  second  emphasised  direct  weight  reduction  by  em- 
ploying numerical-search  methods.  For  comparison  purposes,  all  of  the  candidate 
procedures  were  applied  to  the  inter* edlate-ccnplexity-wing  structure  described 
in  Sjbsection  11.2.1.  The  objective  was  to  resize  an  initial  fully  stressed, 
strength-governed  design  to  achieve  a 3^  increase  in  flutter  speed  with  a minimum 
increase  in  weight.  The  following  paragraphs  summarize  these  study  efforts  and 
provide  the  background  that  led  to  the  final  resizing  algorithm. 

9-2.1  Energy-Based  Optimality-Criteria  Methods 

Although  it  was  recognized  early  in  the  study  that  optimum  resizing  to  in- 
crease flutter  spaed  shotld  aim  toward  achieving  uniformity  among  the  flutter- 
velocity  derivatives  of  all  resized  elements,  it  was  not  evident  that  a simple 
resizing  equation  could  be  formulated  for  satisfying  this  criterion.  Because  of 
the  rather  complicated  nature  of  the  expression  for  these  derivatives,  as  pre- 
viously developed  by  Rudisill  and  Bhatia  in  Reference  9-«l,  it  was  felt  that  it 
would  be  difficult,  if  not  impossible,  to  devise  such  an  equation.  Yet,  it  was 
reasoned  that  for  practical  reuizing  of  very  large  structur.  I idealizations,  some 
flutter-resizing  procedure  that  embodied  the  same  basic  simplicity  as  the  fully- 
stressed-design  approach  for  strength  resizing  should  be  developed,  even  if  it  was 
necessary  to  compromise  the  correct  optimality  criterion.  Along  these  lines, 
several  procedures  based  on  approximate  optimality  criteria  and  simple  energy- 
based  resizing  equations  were  conceived  and  examined.  The  hope  here  was  that 
even  if  the  correct  criterion  was  not  satisfied,  the  resulting  design  would  still 
be  efficient,  although  not  optimum. 

The  two  simplest  versions  of  these  approximate  methods  were  identified  as 
the  "torsi-n  mode  fully  stressed  design"  and  the  "flutter  mode  fully  stressed 
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design"  procedures.  They  were  selected  for  study  mainly  because  they  could  be 
integrated  very  simply  into  an  existing  fully-stresaed-deaign  program  by  the  mere 
addition  of  pseudo-load  conditions,  acting  along  with  the  actual  applied  loadi n& 
conditions,  in  the  resizing  process.  The  first  was  based  on  assumed  proportion- 
ality between  the  fundamental  torsional  frequency  and  the  flutter-speed,  and 
attempted,  in  an  approximate  way,  to  optimize  for  torsional  frequency.  It  used  a 
properly  scaled  set  of  inertial  loads  in  the  torsion  mode  at  a single  additional 
pseudo-loading  condition.  The  second  method  used  a set  of  inerti&l-plus-aero- 
dynaais  loading  conditions  based  on  the  complex  flutter  mode.  In  both  methods, 
resizing  for  strength  and  flutter-speed  requirements  was  performed  simultaneously, 
so  that  interaction  between  the  two  requirements  due  to  internal  load  redistri- 
bution was  achieved  automatically. 

Two  somewhat  more  sophisticated,  but  still  approximate,  procedures  sought 
to  achieve  "uniform  frequency  derivatives  in  the  torsion  mode"  and  "uniform  mean- 
strain-energy  density  in  the  flutter  mode."  The  first  of  these  approaches  still 
relied  on  assumed  proportionality  between  flutter  speed  and  torsional  frequency, 
but  treated  the  frequency-optimization  problem  in  a more  exact  way  than  in  the 
previous  torsion-mode  fully-stressed -design  approach.  The  second  method,  which 
received  previous  attention  by  Siegel  in  Refererue  9*3,  resized  all  flutter-criti- 
cal elements  so  as  to  obtain  equal  values  of  average  strain  energy  per  unit,  weight 
during  a flutter-oscillation  cycle. 

Three  of  the  above  methods  led  to  converged  designs  that  satisfied  the  re- 
quirement of  a 30$  flutter-speed  increase  with  approximately  the  same  increase  in 
structural  weight.  One  method  (the  one  which  attenpted  to  achieve  uniform  fre- 
quency derivatives  in  the  torsion  mode)  behaved  so  poorly  that  resizing  was  aborted 
before  reaching  the  desired  flutter  speed.  To  evaluate  the  final  results  for  the 
converged  designs,  the  flutter-velocity  derivatives  of  the  elements  that  were  re- 
sized to  meet  the  flutter-speed  requirement  were  calculated  and  examined  for  uni- 
formity. The  conclusion  v s discouraging  in  that  no  tendency  toward  uniformity 
existed,  and  no  confidence  could  be  placed  in  any  of  the  methods.  Concurrently, 
results  from  pai-allel  studies  (discussed  in  the  next  subsection)  that  used  numeri- 
cal-search techniques  in  conjunction  with  flutter-velocity  derivatives  confirmed 
the  existence  of  significantly  lighter  designs  than  those  obtained  with  the 
previously  discussed  approximate  optimality-criteria  methods. 


Although  none  of  the  preceding  methods  were  acceptable  in  then selves,  the 
knowledge  gained  from  their  development  became  extremely  valuable.  All  of  the 
methods  employed  energy-type  concepts  — both  in  establishing  their  approximate 
optimality  criteria  and  in  obtaining  simple  and  effective  resizing  formulae  for 
satisfying  these  criteria.  It  then  became  evident  that  by  extending  these  energy 
concepts,  it  war  ncir-r.ble  to  cast  the  expression  for  the  exact  flutter-velocity 
derivatives  of  S.f ■ :,ce  9-2  into  a new  fora  that  was  identifiable  in  terms  of 
generalized  energy  quantities.  Moreover,  it  was  also  possible  to  obtain  a simple 
strain-energy-based  resizing  formula  for  achieving  the  desired  state  of  uniform 
flutter-velocity  derivatives  among  structural  elements. 

When  the  formula  was  applied  to  the  example  wing  structure,  the  final  de- 
sign achieved  the  required  30$  increase  in  flutter  speed  with  a weight  increase  of 
only-about  one-fourth  of  that  required  by  the  previous  methods.  Also,  it  was  ob- 
served that  at  each  intermediate  design  step  in  the  overall  resizing  process,  the 
flutter-velocity  derivatives  of  the  flutter-critical  elements  exhibited  a high 
degree  of  uniformity,  thereby  demonstrating  that  the  resizing  formula  embodied 
excellent  convergence  characteristics.  Nevertheless,  since  the  method  relied  on 
strain-energy-related  quantities,  an  additional  or  more  general  redesign  formula 
was  needed  to  consider  mass-balance  variables.  The  numerical-search  methods  dis- 
cussed next,  as  well  as  the  final  method  selected  for  use  in  FASTOF,  have  the 
capability  for  dealing  with  mass  balance. 

9.2.2  Numerical  Search  Procedures 

Paralleling  the  evaluation  of  the  energy-based  methods  was  the  development 
ar.d  study  of  several  numerical-search  procedures  that  all  employed  the  previously 
referenced  expression  for  the  precise  flutter-velocity  derivatives.  A major  dis- 
tinction between  these  procedures  and  those  already  discussed  is  that  the  numerical- 
search  methols  do  not  rely  on  the  definition  and  enforcement  cf  an  optimality  cri-  ■ 
terion.  Inste*.  i,  concepts  of  travel  in  design  space  are  employed  to  seek  out  a 
near-minimum-we.,-'ht  design  that  satisfies  the  flutter-speed  constraint  without 
compromising  strength  requirements. 

Of  the  four  numerical- search  procedures  studiea,  the  first  two  were  in  the 
category  of  procedures  for  initially  achieving  the  desired  flutter  speed.  The 
latter  two  were  techniques  for  minimizing  structural  weight  after  the  flutter 
speed  target  was  achieved.  A complete  resizing  method  required  that  procedures 
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in  both  categories  be  used  sequentially. 

In  the  first  category,  the  "flutter-velocity-gradient  redesign"  procedure 
was  a straightforward  approach  that  added  increments  of  weight  to  structural  ele- 
ments' in  proportion  to  their  flutter-velocity  derivatives.  The  second  procedure, 
"weight-factor ed-gradient  redesign",  was  a refinement  of  the  first,  wherein  the 
incremental  weight  added  to  each  element  was  factored  by  the  element’s  total 
weight.  This  was  done  in  order  to  arrive  at  the  desired  flutter-speed  constraint 
surface  with  a lower-weight  design  than  thrt  obtained  in  the  previous  method, 
which  tended  to  add  excessive  weight  to  the  lighter  structural  elements.  Goth 
techniques  were  applied  to  the  intermediate-complexity-wing  structure  in  a step- 
by-step  resizing  mode,  arriving  at  the  desired  flutter  speed  with  much  smaller 
increments  of  structural  weight  than  any  of  the  previously  described  approximate 
optimality-criteria  methods. 

To  then  travel  along  the  constraint  surface  to  a minimum-weight  design  point, 
a "biased  tangent"  approach  and  the  "method  of  feasible  directions"  (see  References 
9-4  and  9-5)  were  each  employed  separately,  starting  with  the  last  design  obtained 
by  the  flutter-velocity-gradient  redesign  method.  Both  procedures  led  to  final 
designs  having  essentially  the  same  weight  as  that  achieved  by  the  last  energy- 
based  optimality-criteria  approach  that  aims  for  uniform  flutter-velocity  deriva- 
tives. However,  although  both  procedures  yielded  good  results,  considerable 
difficulty  was  experienced  in  developing  an  efficient  automated  step-size  deter- 
mination procedure. 

In  summarizing  the  findings  of  these  numerical-search  studies,  two  major 
points  should  be  noted.  First,  the  requirement  for  a two-phase  redesign  operation, 
coupled  with  the  problem  of  step-size  determination,  led  to  the  conclusion  that 
these  procedures  are  computationally  inefficient  and  not  readily  amenable  to  com- 
plete automation.  Second,  the  biased-tangent  and  feasible-direction  methods  yield 
the  same  design  a3  that  achieved  by  the  energy-based  resizing  method  that  aims  for 
uniform  flutter- velocity  derivatives,  thus  giving  added  confidence  in  the  superi- 
ority of  this  optimality  criterion. 

9.3  THE  SELECTED  FLUTTER  RESIZIHG  ALGORITHM 

From  the  results  of  the  studies  described  in  the  previous  subsection,  it  was 
concluded  that  the  finally  selected  flutter  resizing  algorithm  should  be  a direct 
rather  than  a two-phase  procedure  that  achieves  a state  of  uniform  flutter-velocity 
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derivatives  for  flutter-critical  elements.  Moreover,  for  the  overall  problem  of 
determining  a near-minimum-weight  design  that  satisfies  both  strength  requirements 
and  a minimum  flutter-speed  constraint  (for  one  critical  flutter  mechanism),  the 
resizing  procedure  should,  aim  toward  a composite  optimality  criterion  with  the 
following  characteristics : 

• Flutter-critical  elements  have  uniform  flutter-velocity 
derivatives  for  the  Mach  number  and  altitude  of  the  prescribed 
critical  flutter  condition, 

• Strength-critical  elements  are  fully  stressed  for  at  least  one 
of  the  specified  design  loading  conditions, 

• Other  elements  are  at  a minimum  or  maximum  manufacturing  gage. 

To  achieve  the  uniform-fit  ;•<, er-velocity-derivative  criterion,  the 
following  resizing  formula  is  used  in  FASTOP: 


where 


old 


new 


oyam^oid 


(^y^iWet 


are  the  ith  design  variable  (structural 
element  or  mass  balance)  weights  before 
and  after  a resizing  step,  respectively, 

is  the  flutter-velocity  derivative  of  the 

XL 

i design  variable  before  a resizing 
step,  and 

is  an  approximation  of  the  desired 
uniform  positive  value  of  the  derivative 
after  a resizing  step,  to  be  discussed 
in  Subsection  9*1*.!. 


This  formula  resulted  from  a simplification  of  the  equation  that  was  used  in  tne 
last  of  the  energy -based  optimality-criteria  methods  discussed  earlier.  It, 
however,  does  not  suffer  from  the  shortcoming  of  its  predecessor,  in  that  it  is 
mathematically  capable  of  dealing  with  mass-balance  as  well  as  structural  design 
variables.  When  this  resiling  technique,  hereafter  referred  to  as  the  "velocity- 
derivative-ratio  method,"  was  applied  to  the  same  example  wing  structure  used  in 
the  previous  study  efforts,  it  displayed  excellent  convergence  characteristics, 
rapidly  leading  to  a design  requiring  the  same  weight  Increase  as  those  obtained 
by  the  best  energy-based  method  and  by  the  numerical-search  procedures.  When  mass- 
balance  design  variables  were  introduced  into  the  example  problem,  the  ability  of 
the  method  to  cope  with  this  type  of  variable,  in  combination  with  structural 
design  variables,  was  also  verified.  The  results  of  further  demonstrations  of 
the  application  of  this  method  in  optimizing  realistic  structural  designs  is 
presented  in  Section  12. 

9.3.1  Calculation  of  Flutter  Velocity  Derivatives  in  FASTOP 

As  stated  previously,  in  the  lat  er  part  of  Subsection  9.2.1,  it  became  evi- 
dent in  studying  the  energy -based  optimality-criteria  methods,  that  it  is  possible  to 
cast  the  analytic  flutter-velocity  derivatives  of  Reference  9-2  into  a new  form  that 
is  identifiable  in  terms  of  generalized  energy  quantities.  The  development  pro- 
ceeds as  follows: 

T 

Consider  the  flutter-mode  vector  (Uj  and  its  associated  row  vector  f VJ  , 
in  the  structures  mathematical  model,  to  be  normalized  such  that 

{V}T([M]  + [A]  ){U)  « 1.  (9.3) 

Rudisill  and  Bhatia's  final  flutter-velocity  derivatives  (Reference  9-2) 
may  then  be  expressed  in  the  following  form: 

gy  y r 

^ - Jji  Re  ({V}T[3[K]/3mi  - aPQ/^HU}) 

la  ({vfCofK]/^  - ajaCHl/SBjCU}) 

M^(V)T(3[A]/3k){U}) 

/ 2 2 
X\y  ^ ([V)T(3CA]/3k)[U})  + 
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(9.4) 


where  [K]  and  [M]  are  the  stiffness  and  mass  iratrices  of  the  structures  mathemat:  cal 
model;  is  the  angular  flutter  frequency;  k.  is  the  reduced  frequency,  bu^/V,, 
in  which  b is  the  reference  wing  semichord;  and  [A]  is  the  complex  aerodynamic 
force  matrix  compatible  with  the- structural  degrees  of  freedom  of  the  surface. 


If  the  ith  element’s  stiffness  and  mass  matrices,  [K^]  and  [Mj],  vary  linearly 
with  it 3 design  variable,  Equation  (9.*0  can  be  manipulated  into  a new  "encrgy- 
densJiy"  form: 


or 


Af 

*1  " 


fVR}T[K1j{UR} 


Af  [[VR}T[MiHUR}  {VI)T[Mi]{UI] 
\2[  Bi  ' ai 

. t tv^w j \ 


3V, 


-■  . SEDf  - KED,  , 


(9-5) 


(9.6) 


where  the  subscripts  R and  I indicate  the  real  and  imaginary  components, 
respectively,  of  the  flutter  vector  and  its  associated  row  rector,  and 
where  C is  a real  quantity  defined  by 


C „ ^ CfVlT(5[A3/dlOfu])  + 2/k  . (9>?) 

Im  ({V)T0[A]/3k)[U}) 

In  Equation  (9.5),  terms  have  been  grouped  into  two  categories.  The 
first,  SEDj,  includes  what  may  be  interpreted  as  a linear  combination  of 
generalized  strain-energy-density  terns.  The  second  category,  KEI^,  contains 
a similar  set  of  generalized  kinetic-energy-density  terms.  This  grouping  of 
terms  is  adopted  in  the  formulation  and  programming  of  the  flutter-velocity 
derivatives  to  etiable  the  FASTOP  user  to  readily  compare  the  separate  contri- 
butions of  each  element’s  stiffneas  and  mass  to  its  derivative.  An  example  of 


10? 


the  usefulness  of  this  feature  occurs  when  identifying  locations  in  the  structure 
where  mass-balance  could  be  effective  in  increasing  a surface's  flutter  speed 
(i.e.,  where  large  negative  KED's  exist).  This  concept  is  discussed  further  in 
Subsection  10.5.2. 

In  order  to  compute  the  derivatives  defined  by  Equation  (9*5) » it  is  necessary 
to  first  transform  the  flutter  vector  and  its  associated  row  vector  from  <:’'e  modal 
coordinates  used  in  the  flutter  analysis  module  to  absolute  structural  coordinates. 
The  different  transformation  relations  required  for  cantilevered  structures  and 
free-free  structures  are  discussed  below. 

(a)  Cantilevered  Structures  - In  cantilever  analyses,  the  transformation 
between  modal  and  structural  coordinates  utilizes  the  normal  node 
shades  £xj  (see  Equation  (7-9))  of  the  dynamics  model,  and  matrix 

which  transforms  displacements  from  dynamics  coordinates  to 
structures  coordinates  (see  Equations  (5*6)  and  (5*7)).  If  lower 
and  upper  case  symbols  are  respectively  used  to  denote  vectors  in 
modal  coordinates  and  structures  coordinates,  the  required  trans- 
formations take  the  form 

{uj  * |q]  {u}  , { V ] T « {vj  T[q]T  (9-8) 

where 

[Q]-[Bf[x]  (9- 9) 

(b)  Free-Free  Structures  - It  may  be  recalled  that  the  basic  vibration 
equation  for  free-free  analyses  (see  Equation  (7.20))  was  cast  in  terms 
of  relative  dynamic  displacements  J^DR  }.  Thus,  {uj  and  {v}  are 
transformed  to  relative  structural  coordinates  by  the  relations 

|vi)-K]i“>  ■ (v™ijT  K]1  • 

where 

h]  - [•»]• 
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Now,  rigid-body  motion  of  structures  points  depends  on  both  the  plug  motion  w 
(see  Equation  (7.19)),-  and  matrix  (Xg]  vhich  defines  those  displacements  for 
unit  plug  notions  (see  Section  7«M*  Accordingly,  |u|  and  |v|  ^ are  transformed 
to  absolute  structural  coordinates  by  the  relations 


Kb.}“[QA]M  *{vab s}T  - T[Qa]T  * 


(9.12) 


W-  M+ MW 


(9.13) 


The  transformation  to  absolute  coordinates  given  by  Equations  (9.12)  and  (9.13) 
is  used  in  the  computation  of  the  kinetic-energy-density  expression  ir  Equation 
(9*5).  Equations  (9.12)  and  (9.13)  can  also  be  used  in  the  computation  of  the 
strain-energy-density  expression.  However,  as  strain-energy-density  is  only  a 
function  of  relative  displacements,  Equations  (9.10)  and  (9.11)  are  used  fcr  that 
computation  since  better  numerical  accuracy  is  achieved. 

9. 3. 1.1  A Theoretical  Consideration.  In  the  computation  of  the  flutter- velocity 

derivatives,  the  normalisation  step  defined  in  Equation  (9.3)  and  the  computation 
of  the  coefficient  C defined  by  Equation  (9.7)  are  actually  carried  out  directly 
in  modal  coordinates  using  [5  J , the  generalised  air  force  matrix  (instead  of  [a]), 
and  the  generalised  mass  matrix  (instead  of  the  structures -model  mass  matrix). 
Moreover,  the  flutter  vector  and  its  associated  row  vector  are  also  initially 
computed  in  modal  coordinates  before  transforming  to  structural  coordinates. 

Now,  in  modal  coordinates,  the  flutter  equation  takes  the  form 


(M-  4 ([».]*[-«])>} 


(9.11*) 


where  [Km]’  [Mm]  *nd  [ft]  tee  generalised  stiffness,  mass  and  earcdyiWc 
matrices,  respectively.  That  is, 

Kb  - [k][q]  , (9.15a) 

Mm  * Qf  [M  * » (9.15b) 

[4Wm 


•Lb,  mA 


where  |k]  , |m]  and  [a|  are  the  structural  stiffness,  mass,  and  aerodynamic 
matrices,  respectively,  and  (q|  is  the  transfonaation  matrix  defined  in  Equation 
(9.9);  for  a free-free  analysis,  matrix  |q^|  (defined  in  Equation  (9-13))  replaces 
[q|  in  Equations  {9. 15).  If  Equation  (9*1*0  Is  used  as  the  starting  point,  the 
derivation  of  the  flutter-velocity  derivative  expression  results  In  an  equation 
identical  to  Equation  (9.4)  except  that  all  vectors  and  matrices  are  in  modal 
rather  than  physical  coordinates;  thus  terms  of  the  form  MT  4Km1  M and 


A»i 

jif  )u}  arise.  In  order  to  evaluate  these  terms,  Equations  (9.15a)  and 
3«i 

(9.15b)  must  be  used.  Now,  the  columns  of  [q]  are  simply  the  vibration  mode 
shapes  expressed  in  structural  coordinates.  If  these  mode  shapes  are  assumed  to 
remain  constant  during  a variation  in  , it  follows  from  Equations  (9.15a)  and 
(9.15b)  that 


(9.16) 


It  is  noted  that  the  terns  on  the  right-hand  side  of  equations  (9*16)  are 
precisely  the  terns  obtained  by  substituting  the  transformations  of  Equations 
(9.8)  or  (9.12)  into  Equation  (9.5)>  as  described  previously.  Thus  it  follows  thut 
the  computation  of  flutter-velocity  derivatives  within  FASTCP  assumes  the  |q] 
matrix  to  be  constant,  i.e.»the  vibration  mode  shapes,  expressed  in  structural 
coordinates,  are  retained  as  the  generalized  displacement  vectors  during  a 
variation  in  the  design  variables. 


9.4  IMPLEMENTATION  OF  THE  COMBINED  STRENGTH  AND  FLUTTER  RESIZING  PROCEDURE 
9.4.1  Determination  of  a Design  Change  for  a Desired  Flutter-Speed  Increment 

In  applying  the  flutter  resizing  formula  of  Equation  (9*2),  an  iterative 
procedure  is  used  to  determine  the  value  of  the  target  derivative,  (BV^/Bm^) target* 
The  procedure  makes  use  of  the  assumption  that  for  small  design  changes,  the 
flutter-velocity  derivatives  may  be  used  to  predict  a change  in  flutter  speed, 

AVf  , as  a linear  combination  of  element  weight  changes,  Am^;  that  is 
pred  * 
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To  start  the  Iteration,  a target  derivative  of  80^  of  the  average  of  all 
positive  derivatives  is  assumed.  Equation  (9*2)  is  then  used  to  determine  a 

new  weight,  m.  , for  each  variable.  The  new  weights  are  then  adjusted  so  that 

new 

no  design  variable  (a)  falls  below  that  required  for  strength,  (b)  violates  mini- 
mum and  maximum  gage  constraints,  or  (c)  is  reduced  by  more  than  the  percentage 
specified  by  user-supplied  "max. -cut"  parameters.  (These  considerations  we  dis- 
cussed more  fully  in  Subsection  9-**.2).  The  resulting  incremental  weights, 

Am ,«  m,  - m,  , are  then  used  in  the  linear  relationship  of  Equation  (9*17) 

1 new  lold 

to  compute  a predicted  velocity  increaent,  dV*  , which  is  then  compared  with 

pred 

the  desired  velocity  increment,  . Two  convergence  criteria  are  used  for 

* 


this  comparison'  satisfactory  convergence  occurs  when 


or  when 


(9-l8a) 


(9.18b) 


where  a and  *2  are  specified  by  the  user.  If  neither  criterion  is  satisfied,  the 
target  derivative  is  automatically  adjusted  within  the  program  (see  discussion 
below)  and  the  entire  procedure  is  repeated  to  obtain  another  trial  redesign.  This 
iterative  process  continues  until  either  a)  at  least  one  criterion  is  satisfied, 
or  b)  it  is  established  that  neither  criterion  can  ever  be  met.  For  the  latter 
case,  wherein  the  design  does  not  change  from  one  trial  resizing  to  the  next,  t.ie 
last  trial  design  is  accepted;  this  situation  may  arise,  for  example,  due  to  the 
presence  of  maximum  gage  constraints. 

The  automatic  adjustment  of  the  target  derivative  from  one  trial  redesign 
to  the  next  proceeds  in  the  following  manner.  Suppose  that  the  first  trial  redesign 
indicates  that  the  predicted  velocity  increment  is  larger  than  the  desired  increment. 


tfo  long  as  this  condition  holds  true  (^V  > ^V.  ),  the  target  derivative  is 

rpred  rdec 

repeatedly  increased  by  10$  of  its  previous  value  in  each  of  the  ensuing  trial 

resizings.  Eventually,  when  the  predicted  increment  becomes  smaller  than  the 

desired  increment  ( ^V’  <•  &V.  ),  the  target  derivative  is  repeatedly 

pred  rdes 

decreased  by  lj  of  its  previous  value ; when  ^V,  again  becomes  larger  than 

pred 

aV_.  . . the  target  derivative  is  successively  increased  by  0.1$  of  its  previous 

-des’  

value.  That  is,  each  time  "crosses"  the  program  changes  both  the 

rpred  aes 

sense  and  magnitude  of  the  target  derivative  adjustments  (+10$,  -1$,  +0.1$,  etc.). 

As  mentioned  above,  the  process  terminates  when  Equation  (o.l8a)  or  Equation  (9.18b) 
is  satisfied,  or  when  the  design  does  not  change  from  one  triad  resizing  to  the  next. 

Clearly,  the  initial  supposition  ( &V.  > ) in  the  first  trial  redesign 

pred  des 

is  immaterial;  if  the  opposite  were  true,  the  starting  target  derivative  would  be 
decreased,  not  increased,  by  1C$. 

9.4.2  Definition  of  "Max. -Cut"  Parameter.*  and  Representation  of  Strength  and 
Manufacturing  Constraints  when  Resizing  for  Flutter 

When  using  the  flutter  resizing  equation,  Equation  (9-2),  the  question  of 
how  to  resize  elements  with  very  small,  or  ever,  negative,  derivatives  must  be 
addressed.  In  such  cases,  it  might  appear  desirable  to  reduce  the  element's  size 
to  its  value  dictated  by  strength  or  minimum-g*ge  requirements.  However,  in  some 
cases  it  has  been  found  that  the  stability  of  the  resizing  procedure  is  improved 
if  the  reduction  in  a structural  element's  size,  in  a single  redesign  step,  is 
constrained  to  a specified  precentage  of  its  previous  value.  Since  it  has  not 
been  found  necessary  to  apply  the  same  restriction  to  mass-balance  design  variables, 
the  user  may  specify  a separate  reduction  factor  fqr  mass  balance  (normally  zero, 
i.e.,  no  restriction).  These  reduction  factors,  referred  to  as  the  "max. -cut"  para- 
meters, are  discussed  further  in  Subsection  10.5.4. 

In  addition  to  these  "max. -cut"  constraints,  strength  requirements  and  mini- 
mum and  maximum  manufacturing  gage  limitations  must  be  considered  when  resizing 
for  flutter.  The  interaction  of  a flutter-speed  constraint  with  strength  require- 
ments is  accounted  for  by  successively  optimizing  for  strength  and  flutter,  with 
the  strength-designed  members  from  the  fully- stressed-design  procedure  (ASOP),  (see 

Section  4)  being  considered  as  minimum  gages  in  the  next  flutter  optimization,  and 
vice  versa. 
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are 


An ' element  * s minimum  and  maximum  manufacturing  gages,  t and  tmpy, 

supplied  by  the  user  during  the  initial  pass  through  the  strength  analysis  and 
redesign  module  of  FASTOP.  The  output  of  that  module  includes  element  gage  size, 
t,  and  stress  ratio,  s (actual  stress/allowable  stress),  which,  along  with  t ^ 

and  t , are  transferred  to  the  flutter  resizing  module.  The  stress  ratio  will 

st&x 

be  very  close  to  unity  for  a stress-critical  element,  but  it  will  be  less  than 

unity  for  a flutter-critical  element,  i.e.,  for  an  element  that  has  been  resized 

for  flutter  (without  encountering  strength  or  minimum  manufacturing  constraints) 

in  a previous  pass  through  the  flutter-resizing  module . The  minimum  strength- 

adequate  gage  size  for  any  element  is  therefore  the  product  of  its  stress  ratio 

and  associated  gage  size,  i.e,,  t • i x t.  Thu3,  when  resizing  for  flutter,  an 

element  is  not  permitted  to  be  sized  below  t , or  t , whichever  is  larger.  Also, 

min  s 

the  element  cannot  be  sized  to  a value  greater  than  t . 


9.^.3  Multiple  Flutter-Redesign  Steps 

FASTOP  allows  the  user  the  option  to  perform  successive  .'lutter-redesign 
steps  without  computing  new  normal  modes  of  vibration  after  each  step.  In  this 
"coupled-mode"  approach,  the  last  set  ci  computed  normal  mode  shapes  are  re- 
tained as  assumed  modes  and  the  changes  in  the  modal  stiffness  and  mass  matrices, 
H “d  K]  , are  given  by  the  following  expressions. 

Case  (a).  Cantilevered  Structure 


(9.19) 


Case  (b).  Free-Free  Structure 


[O  * [yTWU  ■ kl-  [vTWM 


(9.20) 
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Here , and  JjWg ] ere  the  cumulative  incremental  changes  in  the  structures- 
model  atiffhese  and  mass  matrices  since  the  last  normal-mode  computation,  and  Jq]  , 
[^]  and  [QAj  are  toe- transformations  defined  in  Equations  (9*9) » (9*11)  and  (9.13). 

The  adjusted  generalized  stiffness  and  mss  matrices,  I K I andjM  1 , are  then 

- *»  “•'new  <-  “-'new 

obtained  by  adding  the  above  incremental  matrices  to  the  orthogonal  ones  corre- 
sponding to  the  last  set  of  normal  modes;  that  is 


[<L  ■ H ' W 
EvL  ■ M * W ■ 


(9.21) 


(9-22) 


Note  that  those  new  matrices  are  no  longer  diagonal  (i.e.,  the  original  normal- 
mode coordinates  become  "coupled"). 

The  following  section  describes  this  interactive  resizing  process  frem  the 
user's  point  of  view. 


Section  10 


USE  OF  FASTOP  IOR  INTEGRATED  ANALYSIS  AND  DESIGN 

l 

' 10.1  SUMMARY 

FASTOP  is  an  analysis  and  rede.sign  tool  that  can  be  used  to  generate  near- 

• minimus-weight  designs  for  aircraft  structures  subject  to  combined  strength  and 

I flutter-speed- requirements.  Two  basic  programs  are  involved;  one  Is  primarily 

» concerned  with  the  static  strength  problem,  and  the  other  addresses  the  flutter 

condition.  In  the  typical  redesign  process,  each  of  the  programs  is  executed 
; several  times.  Communication  between  the  programs  ensures  that  the  design  re- 

' quirenients  for  one  type  of  constraint  are  not  intentionally  violated  by  the  other. 

| The  material  contained  in  this  section  is  intended  to  guide  a user  through  the 

j entire  redesign  procedure.  Numerous  suggestions  are  put  forward  based  upon  the 

| experience  gained  in  the  solution  of  the  demonstration  problems  to  be  discussed 

in  Section  12. 

* 10.2  ORGANIZATION  CF  FASTOP 
< 

i 

! The  Flutter  And  Strength  Optimization  Package,  FASTOP,  is  comprised  of  a 

! Strength  Optimization  Program,  SOP,  and  a Flutter  Optimization  Program,  FOP. 

! As  shown  in  Figure  10.1,  each  of  these  major  programs  is  organized  on  a modular 

! basis.  The  modules  are  defined  as  follows: 

i 

i SOP  Modules 

A LAM  - Automated  load  Analysis  Module 
; ATAM  - Automated  Transformation  Analysis  Module 

A SAM  - ^Automated  Strength  Analysis  Module 
j A30M  - Automated  Strength  Optimization  Module 

I 

' FOP  Madules 

( 1 - 
j 

i AVAM  - Automated  Vibration  Analysis  Module 

t 

AFAM  - Automated  flutter  Analysis  Module 
AFOM  - Catenated  Flutter  Optimization  Module 
These  acronyms  will  be  used  to  facilitate  the  discussion  throughout  this  section. 
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Figure  10.1  Modular  Organization  of  FACTO P 

10.3  SEQUENTIAL  USE  OF  THE  STRENGTH  AND  FLUTTER  OPTIMIZATION  PROGRAMS 

SOP  and  FOP  may  be  executed  Individually  or,  if  desired,  they  may  be  run 
back-to-back,  as  a multiple  step  job.  Theoretically,  there  is  no  restriction  on 
the  number  of  SOP  and  FOP  steps  that  may  be  executed  in  a given  job.  However, 
practical  considerations  such  as  running  time  and  the  normal  desire  of  the  user 
to  examine  the  output  of  one  program  before  executing  tl.s  next  will  limit  the 
number  of  steps.  It  is  recommended,  therefore,  that  generally  no  more  than  two 
steps  (either  SOP-FOP  or  FOP-SOP)  be  executed  in  a multiple-step  job. 

Figure  10.2  illustrates  the  first  four  steps  in  a typical  redesign  pro- 
cedure. Notice  that  SOP  must  be  the  first  program  executed  in  the  entire  pro- 
cedure and  that  the  two  programs  alternate  thereafter.  As  shown  in  the  figure. 


Start 


Taiie  A - Tine  A — 

SOP 

_ Tape  C 

SOP 

Taoe  C’_ 

Tape  B _ 

FOP 

Tape  B'  __ 

FOP 

Tape  D' 

Tape  D 

figure  10.2  Basic  I/O  Tapes  in  FASTOP 
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data  is  transferred  U ween  programs  by  means  of  magnetic  tapes.  Each  program 
generates  two  basic  output  tapes;  one  is  intended  for  the  next  SOP  step  and  the 
other  is  intended  for  the  next  FOP  step.  The  general  contents  of  these  tapes 
are  briefly  discussed  below. 

Tape  A - This  tape,  referred  to  as  the  "SOTOSO"  tape,  is  generated  by  SOP 
for  use  in  the  subsequent  SOP  step.  Essential  data  such  as  structures  model 
geometry,  boundary  conditions,  applied  loads,  etc.  are  passed  via  this  tape. 

Tape  B - This  "SOTOFO"  tape  ij  generated  by  SOP  and  passed  to  the  next 
POP  step.  Latest  member  sizes  and  the  associated  flexibility  (or  stiffness) 
matrix  are  among  the  data  contained  here. 


Tape  C - POP  accepts  member  sizes  from  the  above  "SOTOFO”  tape,  adjusts 
them  by  the  flutter-resizing  algorithm  and  passes  the  updated  manber  data  to 
the  next  SOP  step  via  this  "FOTOSO"  tape. 


Tape  D - Data  peculiar  to  POP,  such  as  mass  data,  mass  balance  locations, 
etc.,  are  passed  from  one  POP  step  to  the  next  by  means  of  this  "POTOFO"  tape. 

10.4  USE  OF  THE  STRENGTH  OPTIMIZATION  PROGRAM  (SOP) 

At  any  stage  of  the  redesign  procedure,  SOP  can  be  used  in  one  of  three 
possible  modes: 


(a) 

(b) 

(o) 


to  simply  compute  the  dynamics  model  flexibility  matrix  A , or  the 


structural  stiffness  matrix 


H 


to  stress  analyze  the  structure  and  then  compute  |aJ  or  jlCsJ 
to  perform  one  or  more  PSD  (fully  stressed  design)  cycles  and  then 


compute 


[a]  or  Ks 


Use  of  one  of  these  modes  in  any  SOP  step  does  not  preclude  the  use  of  a different 
mode  in  a later  step. 


10.4,1  Conventional  Use  of  SOP 

The  suggested  procedure  is  to  use  mode  (c)  for  the  initial  SOP  step  and  then, 
depending  on  the  degree  of  screr.gth/flutter  interaction  encountered,  to  use  either 
mode  (a),  (b)  or  (c)  in  subsequent  SOP  steps.  Whenever  mode  (c)  is  employed,  the 
user  must  specify  the  number  of  FSD  cycles  to  be  performed.  For  the  initial  SCP 
step,  experience  has  indicated  that  four  redesign  cycles  are  normally  adequate  to 
transform  a preliminary  design  (e.g. , one  having  uniform  gages)  into  a converged 


117 


P 


I 


r 


t 


fully  stressed  design.  However,  a single  cycle  should  be  sufficient  in  later  SOP 
steps.  In  fact,  the  results  of  the  demonstration  cases  described  in  Section  1 2 
indicate  that  strength  redesign  could  have  been  bypassed  in  all  intermediate  re- 
design cycles  (mode  (a))  with  a final  strength  resizing  near  the  end  of  the 
flutter  redesign  process.  The  user  must  also  specify,  in  the  initial  SOP  step, 
which,  if  any,  of  the  total  set  of  elements  are  to  be  permanently  excluded  from 
the  strength  resizing  process.  This  elimination  becomes  necessary,  for  example, 
when  some  elements  are  adequately  modeled  to  simulate  the  stiffness  of  a structure 
'but  \re  inadequately  modeled  for  stress  analysis  and  redesign  - as  when  honeycomb 
core  is  modeled  with  rib  and  spar  webs.  Actual  exclusion  of  an  element  is  effec- 
ted by  simply  setting  its  minimum  and  maximum  allowable  gages  equal  to  its  initial 
(desired)  gage.  Note  that  since  these  same  allowable  gages  are  transmitted  to 
POP,  these  elements  are  in  fact  withdrawn  fran  the  entire  redesign  procedure. 

When  mode  (b)  is  used,  strength  resizing  does  not  take  place,  but  the  ele- 
ment stress  ratios  (actual  stress/allowable  stress)  are  updated  to  define  revised 
strength  gage  requirements  of  flutter  design  variables  for  use  in  the  subsequent 
POP  step.  This  enables  POP  to  avoid  violating  strer.gth  requirements  when  resizing 
for  flutter. 

10.4.2  Use  of  SOP  with  Strength-Governed  Designs  not  Generated  by  SOP 

One  other  important  use  of  FASTOP  occurs  when  SOP  is  not  intended  to  be 
used  for  stress  analysis  or  redesign.  This  condition  arises  when  u strength-ade- 
quate design,  generated  using  criteria  other  than  those  embodied  in  SOP.  is  found 
to  be  flutter  deficient.  Under  these  circumstances,  mode  (a)  should  be  uted  in 
all  SOP  steps  ar.d,  in  the  first  SOP  step,  the  minimum  allowable  gage  of  each  ele- 
ment should  be  set  equal  to  its  initial  gage.  This  will  ensure  that  any  subse- 
quent flutter  redesign  in  the  POP  steps  will  lead  to  a flutter  adequate  design 
in  which  member  gages  are  nowhere  less  than  those  of  the  Initial  design. 

10.4.3  Initial  Design 

Although  resizing  may  have  occurred  in  the  first  SOP  step,  the  resulting  de- 
sign of  this  step  - not  the  input  preliminary  design  - is  referred  to  as  the 

initial  design.  It  is  this  initial  strength-adequate  design  which  may  be  flutter 

» 

deficient  and  require  subsequent  flutter/strength  resizing  to  achieve  an  adequate 
design.  If  so,  all  changes  in  design  parameters,  such  as  flutter  speed  ar.d  total 
weight  are  defined  with  respect  to  the  initial  design. 


10.5  USE  OF  THE  FLUTTER  OPTIMIZATION  PROGRAM  (FOP) 

There  are  five  possible  nodes  in  which  FOP  can  be  used: 

(a)  to  perform  a vibration  analysis 

(b)  to  perform  a flutter  analysis 

(c)  to  perform  both  vibration  and  flutter  analyses 

(d)  to  perform  vibration  and  flutter  analyses  and  also  compute 
flutter  velocity  derivatives 

(e)  to  perform  vibration  and  flutter  analyses,  compute 
flutter  velocity  derivatives  and  then  perform  one  or  more 
-cycles  of  flutter  redesign  (each  redesign  cycle  is  followed  by 
a "coupled  mode"  flutter  analysis  and  computation  of  flutter 
velocity  derivatives) 

It  is  not  unlikely  that  each  of  these  modes  will  be  exercised  at  some  point  in 
the  entire  redesign  procedure. 

10.5.1  Determination  of  Critical  Flight  Condition 

After  the  first  SOP  step,  the  user  must  determine  if  the  initial  design 
is  flutter  adequate.  The  first  requirement  is  to  obtain  a realistic  dynamic  mass 
matrix  for  the  initial  design.  This  task  is  normally  the  responsibility  of  a 
weights  engineer  who  must  account  for  the  weight  of  detailed  structural  items 
(rivets,  fittings,  etc.),  and  other  items  (e.g.,  engines,  fuel,  actuators,  ex- 
ternal stores),  as  well  as  the  idealized  structural  weight.  POP,  however,  does 
have  the  additional  capability  of  automatically  generating  a dynamic  mass  matrix 
by  using  the  idealized  structural  weights  in  conjunction  with  nonoptimum  factors 
(see  Section  10. 5. 3)  and  additional -mass  data  supplied  by  the  user.  The  next  step 
is  to  perform  a series  of  vibration  and  flutter  analyses  for  various  flight  condi- 
tions in  order  to  determine  the  critical  condition  from  the  standpoint  of  flutter. 
The  vibration  analysis  need  only  be  done  or.ce  for  each  weight  condition,  either 
independently  (mode(a)),  or  in  conjunction  with  the  first  flutter  analysis  for 
that  weight  condition  (mode  (c)).  In  any  event,  the  vibration  data  is  saved,  on 
tape  and  used  in  subsequent  flutter  analyses  (mode  (b)). 
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10.5-2  Initial  Mass  Balance  Data 


If  a design  Is  flutter  deficient  and  the  user  Intends  to  Include  mass 
balance  In  a redesign  study,  the  location  and  an  initial  value  of  each  mass  bal- 
ance oust  be  specified.  Experience  has  shown  that  the  choice  of  this  initial 
Cf.ta  can  have  a significant  effect  on  the  subsequent  redesign  behavior.  For 
example,  a case  was  encountered  in  which  the  effectiveness  of  a mass  balance  de- 
pended strongly  upon  Its  initial  value;  small  initial  values  were  not  useful  (nega- 
tive flutter-velocity1  derivatives)  and  were  subsequently  eliminated  from  the  de- 
sign, while  larger  values  were  found  to  be  very  effective.  Accordingly,  it  may  be 
profitable  to  do  a small  separate  mass-balance  study  (mode(d)),  varying  the  initial 
values  and  locations.  The  resulting  flutter-speed  increments  and  mass-balance 
velocity  derivatives  should  guide  the  urer  in  the  selection  of  good  initial  mass- 
balance  data.  Note  that  the  velocity  derivatives  of  the  structural  elements  also 
contain  useful  information.  Specifically,  the  distribution  of  the  hinetlc-enercy- 
denslty  (KED)  components  of  these  derivatives  is  a direct  measure  of  how  flutter 
speed  will  be  affected  by  small  mass  increments  throughout  the  structure; 
mass  increments  are  most  beneficial  in  regions  of  large  negative  KED  components. 

It  is  also  suggested  that  the  initial  mass-balance  data  include  a number  of 
selected  "dummy”  locations  where  zero  values  of  mass  balance  are  specified.  By 
means  of  this  contrivance,  the  user  essentially  "reserves  the  right"  to  introduce 
real  values  of  mass  balance  at  these  locations  in  any  subsequent  FOP  step;  that  is, 
after  the  initial  FOP  step,  the  user  will  have  the  option  of  changing  the  mass 
balance  value  at  any  location  specified  in  the  initial  data. 

Finally,  one  last  point  must  be  made  regarding  mass-balance  locations.  Con- 
sider an  initial  mass  balance,  m-j.  Unless  the  automatic  mass  generator  option  is 
being  employed,  the  user  must  insert  mj  directly  into  the  initial  dynamic  mass 
matrix  at  the  translational  degrees  of  freedom  associated  with  its  dynamic  model 
node  point,  "d”.  However,  since  all  redesign,  structural  as  well  as  mass  balance, 
is  accomplished  in  the  strictures  model,  any  Incremental  mass  balance  Amj  is 
first  assigned  to  structural  node  point  "a"  (specified  by  the  user)  and  then  trans- 
formed to  the  dynamics  model  by  means  of  the  transformation  defined  in  Equation 
(5.10).  It  is  essential  that  the  transformation  of  abj  from  node  s be  made 
directly  to  node  d and  to  no  other  dynamic  node.  The  user  should  Keep  this  re- 
quirement in  mind  when  creating  the  force  beaming  table  that  prescribes  the  beaming 
of  unit  loads  from  dynamic  nodes  to  structural  nodes.  This  problem  does  not  arise 


120 


when  the  automatic  mast  generator  option  is  employed  because  both  the  initial  and 
incremental  mass  balance  are  then  specified  in  the  structures  model  and  trans- 
formed to  the  dynamics  model. 

10.5.3  Flutter  Redesign  Slcmer.ts/.’fonoptimua  Factors 

After  the  initial  mass-balance  data  has  been  specified,  the  user  should 
indicate  which,  if  any,  of  the  total  set  of  structural  elements  are  to  be  perma- 
nently excluded  from  the  flutter  resizing  process.  As  a minimum,  those  elements 
which  were  previously  eliminated  free  resizing  in  SOP  should  likewise  be  ex- 
cluded in  FOP.  The  user  may  now  also  introduce  nonoptimum  factors  for  use  in 
the  SOF/POP  redesign  process  - unless,  of  course,  these  factors  were  already 
introduced  in  conjunction  with  the  automatic  mass-generator  option,  ifonoptimrun 
factors  are  intended  to  account  for  the  fact  that  there  is  incremental  weight, 
other  than  that  of  the  primary  structure,  associated  with  the  redesign  of  each 
element;  for  example,  when  redesigning  an  element  with  a nor. optimum  factor  of 
1.2,  the  true  incremental  weight  would  be  taken  to  be  twenty  percent  larger  than 
the  computed  incremental  structural  weight  of  the  finite  element. 

10. 5. 4 Basic  Parameters  for  Automated  Flutter  Redesign 

Whe-  ever  mode  (e)  is  exercised  in  a FOP  application,  the  user  must  specify 
input  data  to  control  the  number  of  redesigns  to  be  accomolished  in  the  step, 
the  desired  flutter  speed  step  sizes,  etc.  Same  important  parameters  are  dis- 
cussed below. 

10. 5.4.1  Flutter  Band.  Let  Vf  and  "v^jej  denote  the  current  flutter  speed  and 

the  desired  final  flutter  speed,  respectively.  Whereas  speeds  much  larger  than 
Vfdes  are  undesirable  because  of  the  weight  penalty  associated  with  the  extra 
speed,  values  nominally  in  excess  of  -re  considered  to  be  acceptable. 

Thus,  the  user  must  specify  both  Vfdes  an>'  an  additional  parameter  in  order 
to  define  a "band"  of  acceptable  flutter  speeds;  it  is  suggested  that  this  band- 
width not  exceed  one  percent  of  Vfdes,  in  which  case  the  band  will  extend  from 

vfdes  t0  1*01  'fdes* 

10.5.4.2  Step  Size/ 'formal  vs.  Coupled  Modes.  A step  size  parameter,  :3AR,  and 
a parameter  defining  the  maximum  permissible  number  of  automatic  redesioi  steps, 
I.TIX,  are  also  required  by  FOP.  The  program  first  undertakes  to  raise  (or  lower) 
the  flutter  speed  from  V to  a value  in  the  center  of  the  flutter  band,  V*,  in 
:i3AR  approximately  equal  speed  increments;  thereafter,  each  successive  redesign 

121 


JU* 


attempts  to  get  to  V*  directly.  Execution  continues  in  this  manner  until  the 
design  converges,  or  until  NFIX  redesigns  have  been  effected.  The  first  flutter 
analysis  in  the  foregoing  procedure  utilizes  normal  modes  of  vibration,  but  all 
subsequent  analyses  for  modified  designs  are  based  on  a coupled  mode  approach 
utilizing  the  original  normal  modes  together  with  modified  generalized  mass  and 
stiffness  matrices  including  off-diagonal  terms.  Experience  has  indicated  that 
coupled  mode  results  are  somewhat  unreliable,  and  it  is  therefore  recommended 
that  generally  only  one  redesign  step  be  performed  in  each  FOP  application,  i.e., 
the  user  '.s  advised  to  set  NFIX  * 1.  This  restriction  is  not  very  severe  because 
studies  have  shown  that  good  results  can  be  obtained  with  fairly  large  step  sizes 
and  consequently  a small  number  of  redesign  cycles,  when  the  flutter  analyses  and 
velocity  derivative  computations  are  based  on  revised  normal  modes.  For  example, 
the  flutter  speed  of  the  Intermediate-Complexity  Wing  (see  Section  12)  was  in- 
creased by  thirty  percent  to  a nearly  converged  optimum  design  in  just  four  single- 
redesign  FOP  steps. 

10. 5. ^.3  'Max. -Cut"  Parameter.  In  preliminary  studies  of  the  aforementioned 
Intermediate-Complexity  Wing,  it  was  found  that  a few  elements  were  undergoing 
severe  fluctuations  in  gage  size  from  one  flutter  redesign  to  the  next.  The 
phenomenon  was  attributed  to  the  fact  that,  due  to  the  coarseness  of  the  model, 
the  load  pi.ths  were  very  sensitive  to  design  changes.  This  stability  problem 
was  resolved  by  simply  not  allowing  any  gage  size  to  be  reduced  by  more  than 
twenty- five  percent  in  any  single  redesign,  i.e.,  the  "max-cut"  parameter,  D, 
was  set  equal  to  0.75.  Difficulties  of  this  sort  did  not  occur  in  the  redesign 
study  of  the  all -movable  stabilizer  - for  which  a very  detailed  model  existed. 
Indeed,  the  entire  resizing  procedure  progressed  very  smoothly  and  no  restriction 
had  to  be  imposed  on  gage  size  reduction  (D  * 0.0).  Accordingly,  the  user  is 
advised  to  begin  the  redesign  procedure  with  D » 0.0;  then,  if  gage-size  insta- 
bilities appear,  the  "max-cut"  parameter  can  be  adjusted.  Note  that  the  user 
must  specify  a separata  "max-cut"  parameter,  DEAL,  for  mass  balance  variables. 

10.5.5  Termination  of  Redesign  Process 

It  has  already  been  pointed  out  that  an  efficient  design  (close  to  the 
optimum)  can  usually  be  achieved  with  fairly  large  step  sizes  provided  the  single- 
ctep,  normal-mode  approach  is  employed.  Once  a flutter  adequate  design  has  been 
obtained,  the  user  should  avoid  excessive  iterations  within  the  flutter  band 
while  striving  to  effect  a condition  of  uniform  flutter-velocity  derivatives. 
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Experience  has  shown  that  such  an  approach  can  expend  large  amounts  of  ccnputing 
time  for  small  improvements  in  design  efficiency.  It  is  recommended,  therefore, 
that  the  redesign  process  he  terminated  when  weight  reductions  from  two  or  three 
successive  iterations  are  no  longer  significant  from  an  engineering  point  of  view. 


123 


Section  11 


i 


EXAMPLES  n F ANALYSIS  RESULTS 

11.1  SUMMARY 

Calculation#  to  demonstrate  the  capability  and  versatility  of  the  FASTOP 
analysis  package  have  been  performed  for  three  typical  designs  - a prelimin- 
ary structures  model  of  a simple  cantilevered  wir.g  and  detailed  models  of  an 
all-movable  stabiliser  and  a wing  with  a pylon-mounted  external  store. 

The  structures  of  all  three  surfaces  were  represented  by  finite-element 
idealizations  and  the  degrees  of  i.eedom  of  the  structures  models  were  re- 
duced to  a lesser  number  of  degrees  of  freedom  required  for  the  dynamics 
models  by  e force  beaming  transformation  procedure.  The  types  of  calcula- 
tions performed  include  subsonic  and  supersonic  aerodynamic  load  predictions. 
Inertial  load  predictions,  load  beaming  to  structural  node  points,  determina- 
tion of  internal  member  loads  and  stress  ratios,  formation  of  dynemics  model 
flexibility  matrices,  vibration  mode  analyses,  and  subsonic  and  supersonic 
flutter  analyses  using  both  the  k and  p-k  solution  procedures.  Typical  re- 
sults obtained  for  the  three  demonstration  problems  are  presented  and  dis- 
cussed. 

11 .2  DESCRIPTION  OF  DEMONSTRATION  PROBLEMS 

11.2.1  Structures  Model  of  the  Intermediate  Complexity  Wing 

The  structural  idealization  for  the  so-called  "intermediate  complexity  wing" 
is  the  simplest  of  the  three  demonstration  cases  selected  and  is  representative 
of  a typical  preliminary  design  configuration-  The  all-aluminum  two-cell  wing 
box,  illustrated  in  Figure  11.1,  is  modeled  using  100  finite  elements. 

Membrane  element.!  are  used  to  represent  the  wing  covers,  shear  panels  repre- 
sent the  spar  and  rib  webs,  and  bar  elements  are  introduced  between  upper  and 
lower  cover  node  points.  The  wing  root  is  built-in  by  fully  constraining  all 
structural  nodes  on  the  root  boundary.  The  structure  has  a total  of  140  degrees 
of  freedom  comprising  3 trar.rletional  degrees  of  freedom  at  each  structural 
node  point. 
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11.2.2  Dynamics  Model  of  the  Intermediate  Complexity  Wing 

The  150  degrees  of  freedom  of  the  structures  model  are  transformed  to 
39  dynamic  degrees  of  freedom  using  the  force  beaming  procedures  described 
in  Section  5.  The  selected  dynamic  degrees  of  freedom,  indicated  in  Figure 
11.2,  include  out-of-plane  (z  direction)  displacements  at  all  dynamic  nodes, 
these  nodes  being  located  between  adjacent  upper  and  lower  cover  structural 
nodes  in  the  mid-plane  of  the  wing.  The  overhanging  panels,  representing 
the  mass  and  moment  of  inertia  of  structure  external  to  the  wing  box  are 
also  included  in  the  dynamics  model.  Their  inertia  effects  are  included 
by  introducing  rotational  degrees  of  freedom  of  each  panel  about  an  axis 
parallel  to  its  attaching  structure  i.e.,  front  beam,  res"  beam,  or  tip  rib. 
The  local  axis  system  for  each  panel  is  indicated  in  Figure  11.2. 

11.2.3  Structures  Model  of  the  All-Movable.  Stabilizer 

The  all-movable  stabilizer  model,  shown  in  Figure  11.3,  is  representa- 
tive of  a complex  detailed  design  configuration.  A total  of  891  finite 
elements  are  used  to  model  the  stabilizer  surface  including  its  pivot  and 
actuator  restraints.  (The  details  of  the  pivot  restraints  and  the  actuator 
are  omitted  from  the  figure  to  preserve  clarity  of  presentation) . It  should 
be  noted  that  the  inner  and  outer  stabilizer-to-pivot  support  points  are 
modeled  in  the  mid-plane  of  the  surface,  enclosed  by  structural  nodes  h8l, 

U83,  ^63,  l»6l  and  structural  nodes  379,  381,  3^9,  3^7  respectively.  The 
stabilizer  construction  consists  of  titanium  covers,  modeled  as  membrane 
elements,  with  an  aluminum  honeycomb  core.  The  honeycomb  core  is  modeled 
3S  spar.wise  and  chordwise  shear  panel  elements  with  stiffness  properties 
representative  of  the  actual  honeycomb  structure.  The  model  has  a total  of 
1172  degrees  of  freedom, 

11,2.1*  Dynamics  Model  of  the  Ill-Movable  Stabilizer 

The  force  beaming  procedure,  described  in  Section  5,  is  used  to  transform 
the  1172  structural  degrees  of  freedom  of  the  stabilizer  to  92  dynamic  decrees 
of  freedom.  Vertical  (out-of-plane)  degrees  of  freedom  are  specified  at  the 
73  dynamic  node  points  shown  in  Figure  11. U,  and  additional  rotational  decrees 
of  freedom  are  specified  at  overhanging  panel  points,  designated  as  points 
1 through  1?  and  65  through  72. 


Figure  11.2  Intermediate  Complexity  Wing  Dynamics  Model 


11.2.5  Structures  Model  of  the  Wing- with- Store 

A second,  but  significantly  different  detailed  design  configuration  is 
the  wing  with  external- store  shown  in  Figure  11. 5.  The  structural  idealiza- 
tion of  the  multi-spar  aluminum  wing  box  uses  warped  quadrilateral  membrane 
elements  to  represent  the  wing  covers  and  quadrilateral  shear  panels  for  rib 
and  spar  webs.  The  semi-wing  is  modeled  from  the  tip  to  the  airplane  center- 
line  and  symmetric  structural  boundary  conditions  are  specified  for  all  nodes 
in  the  plane  of  symmetry.  The  wing-to- fuselage  connection  is  accomplished 
with  vertical  shear  attachments  at  nodes  37  and  UT  and  a drag  attachment  at 
node  1*7.  A wing-store  pylon,  modeled  with  beam  elements,  is  attached  to  the 
wing  at  nodes  175,  176,  and  l8l,  182.  The  structural  model  haa  a total  of 
602  elements  and  885  degrees  of  freedom. 

11.2.6  Dynamics  Model  of  the  Wing-with-Store 

The  dynamics  model,  illustrated  in  Figure  11.6,  is  schematically  similar 
to  the  two  previous  model;?.  Hov.’ever,  for  the  wing-with-store  there  is  an 
additional  requirement  to  include  lore-and-aft  degrees  of  freedom  at  every 
node  to  account  for  dynamic  coupling  between  wing  pitch  and  store  trans- 
lation. The  store  is  allowed  5 degrees  of  freedom  and  the  dynamics  model 
contains  a total  of  136  degrees  of  freedom. 

11.3  DISCUSSION  OF  ANALYSIS  RESULTS  OBTAINED  WITH  FASTOP 

11.3.1  Loads  Analysis 

Data  pertaining  to  the  basic  aerodynamic  characteristics  of  each  surface, 
plus  a summary  of  flight  conditions  for  which  aerodynamic  load  distributions 
were  computed  using  FASTOP,  are  presented  in  Table  11.1.  Two  examples  of 
pressure  distributions  are  presented  for  the  intermedia te-complexi typing 
demonstration  jase.  The  first  example  (Figure  11. 7)  shows  the  subsonic 
pressure  distribution  for  the  Mach  0.9  flight  condition  noted  in  Table  11.1. 
The  pressures  were  computed  by  the  vortex-lattice  aerodynamics  routine  des- 
cribed in  Section  3.  The  pressure  distribution  for  the  Mach  2.0  flight  condi- 
tion, computed  by  the  supersonic  source  distribution  aerodynamics  routine,  is 
presented  in  Figure  11.8.  In  the  supersonic  case,  the  pressure  is  almost 


Figure  11.5  Wing-vith-Store  Structures  Model 
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uniform  in  the  mid-span  region  due  to  the  almost  two-dimensional  flow  char- 
acteristics at  high  Mach  number.  The  abrupt  pressure  reductions  in  the  in- 
board and  tip  regions  of  the  wing  occu.’  at  locations  where  the  root  and 
tip  Mach  lines  intersect  the  chords  along  which  pressures  have  been  computed. 
Symmetric  boundary  conditions  have  been  specified  for  both  cases,  although 
loads  may  be  optionally  calculated  for  antisymmetric  or  asymmetric  conditions. 

FASTOP  was  also  used  to  compute  inertial  load  distributions  for  the  in- 
termediate complexity  wing  for  translational  and  angular  accelerations  which 
could  correspond  to  the  aerodynamic  flight  conditions  of  Table  11.1,  The 
weights  model  used  for  this  purpose  closely  .resembled  the  dynamics  model  pre- 
sented in  Figure  11.2,  except  that  overhanging  panel  mass  and  inertia  prop- 
erties were  specified  at  panel  c.g. 's  in  an  unswept  axis  system.  The  results 
of  the  inertial  load  analysis  are  in  the  form  of  inertial  forces  and  moments 
in  the  weights  grid. 

11.3.2  Structural  Analysis 

The  structural  analysis  module  was  used  for  each  example  structure  to 
determine  its  level  of  strength  adequacy  and  to  establish  its  flexibility 
characteristics  for  subsequent  vibration  analysis.  Finite  element  sizes  in 
each  structural  idealization  were  selected  to  be  reasonable  but  are  regarded 
as  only  preliminary. 

The  previously  computed  aerodynamic  and  inertial  applied  loads  were 
transformed  to  the  structures  models  by  using  the  procedures  of  Section  5, 
and  ratios  of  maximum  working  stress  to  allowable  stress  (stress  ratios) 
were  computed  for  every  finite  element  in  each  model.  Since  these  stress 
ratios  are  based  on  initial  element  sizes  (for  a non- fully- stressed  design), 
they  are  of  little  significance  except  to  demonstrate  the  proper  working  of 
the  program. 
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In  the  same  analysis  module,  the  structural  stiffness  matrices  were 
assembled.  These  were  subsequently  transformed  into  dynamics  model  flexi- 
bility matrices  by  again  utilizing  the  transformation  procedures  of  Section  5. 

11.3.3  Vibration  Analysis 

Vibration  analysis  mode  shape  plots  for  the  three  demonstratici  struc- 
tures are  presented  in  Figures  11. 9-11. 11.  These  plots  were  generated  using 
the  CALCCMP  plotter  and  associated  software  routines  which  are  part  of  the 
FASTOP  system.  In  each  case  '..e  plotted  data  consists  of  modal  displace- 
ments of  the  surface  dynamic  node  points.  A summary  of  the  modal  character- 
istics of  each  demonstration  structure  appears  in  Table  11.2.  Some  parti- 
cular features  of  the  vibration  analysis  results  follow. 

The  inteimediate-complexity-ving  mode  shapes  are  typical  for  a simple 
cantilevered  structure.  Wing  torsion  (Figure  11.9b),  which  is  the  second 
mode,  has  a node  line  running  spenwise,  in  close  proximi*y  to  the  trailing 
edge.  Chordwise  bending  is  evident  in  modes  2-t. 

The  stabilizer  mode  shapes  (Figures  ll.lOa-d)  indicate  the  presence 
of  root  rotational  motion  due  to  the  pivot  and  actuator  support  flexibil- 
ities. Because  of  this  flexibility,  the  stabilizer  has  both  a pitch  mode, 
in  which  the  stabilizer  rotates  about  its  pivot  (Figure  11.10b),  and  a tor- 
sion mode,  in  which  the  surface  twists  with  virtually  no  pivot  rotation 
(Figure  ll.lOd). 

The  mode  shape  plots  i or  the  wing-with- store  example  are  presented  in 
Figures  ll.lla-f.  The  fore-aft  deflections  for  the  wing  leading  edge  are 
included  on  a separate  reference  line  located  in  the  lower  portion  of  esch 
figure.  The  hand  plotted  store  motions,  both  translational  and  rotational, 
are  shown  with  respect  to  the  undisplaced  store  position.  The  initial  re- 
quirement to  make  the  wing  flutter-critical  was  accomplished  by  defining 
relatively  high  store  mass  and  inertia  properties.  The  re  : ting  high 
stoi^  pitch  inertia  causes  significant  wing  torsion  in  the  store  pitch  mode, 
shown  ir  Figure  11.11c. 
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►.ermediate  Complexity  Wing  Vibration  Mode,  Sheet  1 of  >» 


Figure  11.9  Intermediate  Complexity  King  Vibration  Mode,  Sheet  2 of  it 
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Figure  11.11  Wing- with- Store  Vibration  Mode,  Sheet 


Figure  11.11  Wi ng-vith-Store  Vibration 


11. 3. b Flutter  Analysis 
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All  cf  the  flutter  plots  presented  in  this  section  were  obtained  using 
& CALC0I-5P  plotting  routine,  which  is  an  integral  part  of  FASTOP.  It  should 
be  mentioned  that  this  routine  plots  the  data  points  using  symbol  identifi- 
cation, while  the  curve  fairing  is  & hand  process. 

The  intennediate*<omplexityving  flutter  analysis  for  Mach  0.6,  sea 
level,  was  performed  using  both  the  doublet-lattice  and  assumed-pressure  - 
function  procedures.  The  planform  was  idealized  by  72  aerodynamic  panels 
in  the  doublet-lattice  approach  and  the  chordwise  and  spanwise  polynomials 
assumed  for  the  pressure- function  method  were  3rd  and  6th  degree,  respec- 
tively. Flutter  results  for  both  aerodynamic  representations  using  the  k- 
method  of  solution  are  presented  in  Figures  11.12  and  11.13.  An  additional 
analysis  that  uses  the  doublet-lattice  approach  and  the  p-k  method  of  solu- 
tion is  presented  in  Figure  11. lb.  All  three  analyses  yield  almost  identical 
flutter  speeds.  However,  the  subcritical  frequency  and  damping  trends  differ. 
In  particular,  the  frequency  coalescence  of  the  bending  and  torsion  modes, 
which  causes  the  flutter  instability,  is  more  evident  in  the  results  from  the 
p-k  method  of  solution. 

The  flutter  analysis  of  the  all-movable  stabilizer  was  performed  using 
the  Mach-box  method  for  a Mach  number  of  1.6  and  an  altitude  of  30,000  feet. 
The  aerodynamic  surface  was  represented  by  33b  rectangular  boxes  and  the  free- 
stream  diaphragm  region,  defined  by  the  aft  Mach  line  from  the  outboard  tip 
of  the  ieedir.g  edge  and  the  forward  Mach  line  from  the  outboard  tip  of  the 
trailing  edge,  was  represented  by  two  diaphragm  boxes.  The  small  number  of 
required  diaphragm  boxes  concentrated  in  the  tip  region  is  explained  by  the 
fact  that  both  the  leading  and  trailing  edges  of  the  stabilizer  are  supersonic 
at  Mach  1.6.  The  flutter  analysis  results,  presented  in  Figures  11. 15a, b 
show  that  the  flutter  mechanism  is  caused  by  a coupling  between  the  stabilizer 
bending  and  pitch  modes. 

Flutter  analyses  of  the  virg-with-store,  using  the  doublet-lattice 
method,  are  presented  in  Figures  11.16  a,  b.  The  results,  computed  for 
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Mach  0.8,  sea  level,  clearly  show  a number  of  store  modes  with  veiy  low  aero- 
dynamic damping.  The  lowest  flutter  instability  occurs  at  330  knots  equiva- 
lent airspeed  and  is  caused  by  coupling  between  the  first  wing  bending  mode 
and  the  store  pitch  mode.  The  damping  trend  of  the  critical  root  is  char- 
acteristic of  the  "grazing"  instabilities  encountered  in  store  flutter 
problems. 
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b.  Damping  vs  Airspeed 

Figure  11.12  Intermediate  Complexity  Wing  Flutter  Analysis  (Doublet  Lattice 
Aerodynamics,  k Solution),  Sheet  2 of  2 
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Figure  11.13  Intermediate  Complexity  Wing  Flutter  Analysis  (Assumed  Pressure  9 
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Figure  11. 1U  Intermediate  Complexity  Wing  Flutter  Analysis  (Doublet  Lattice 
Aerodynsmics * Solution) t Sheet-  1 of  2 
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b.  Damping  vs  Airspeed 


Figure  11.1b  Intermediate  Complexity  Wing  Flutter  Analysis  (Doublet  Lattice 
Aerodynamics,  p-k  Solution),  Sheet  2 of  2 
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b.  Damping  vs  Airspeed 


Figure  11.16  Wing- with- Store  Flutter  Analysis  (Doublet  Lattice 
Aerodynamics,  k Solution),  Sheet  2 of  2 


Section  12 


l EXAMPLES  OF  OPTIMIZATION  RESULTS 

| 12.1  SUMMARY 

V 

% 


t Calculations  to  demonstrate  the  redesign  capability  of  FASTOP  have  been  per- 

| formed  using  the  intermediate  complexity  wing,  all-movable  stabilizer,  and  wing- 

| vith-store  models.  The  characteristics  of  these  models  are  described  in  Sec- 

; tion  11.  The  csntilever-structure  dynamics  model  of  the  wing- with- store,  des- 

? cribed  in  Section  11,  was  converted  to  a free- free  model  for  flutter  redesign. 

t 

\ This  was  done  to  fully  demonstrate  FASTOP's  redesign  capabilities  for  a con- 

1 figuration  where  dynamic  interaction  between  the  wing  and  fuselage  might  be 

| important.  The  details  of  the  changes  that  were  made  to  the  dynamics  model  are 

j included  in  the  discussion  of  wing-with-store  redesign  results. 

j Initially,  the  strength  optimization  portion  of  FASTOP  was  used  to  obtain 

< a fully  stressed  design  (FSD)  for  each  demonstration  structure.  In  order  to 

i accomplish  this  calculation,  it  was  necessary  to  specify  flight  design  load  condi- 

i 

i tions  for  the  loads  analysis  module,  initial  member  sizes  in  the  strength  optimi- 

zation module,  and  force  beaming  data  required  for  the  transformation  analysis 
module.  The  latter  ^joup  of  data  was  needed  to  create  transformation  matrices 
from  the  aerodynamics  and  dynamics  mathematical  models  to  the  structures  mathe- 
i matical  model.  Thus  in  the  initial  FASTOP  application  the  design  loads  were  com- 

; puted,  these  were  then  beamed  to  the  structures  model  end  the  structure  was  re- 

i sized  for  strength.  The  analysis  terminated  with  computation  of  a flexibility 

i 

1 matrix  for  the  dynamics  model.  The  subsequent  redesign  to  achieve  a specified 

flutter  speed  improvement  was  achieved  through  multiple  sequential  submissions  of 
■ the  two  major  programs  of  FASTOP  - the  structural  optimization  program  and  the 

flutter  optimization  program.  Except  when  otherwise  uoted,  strength  redesign  was 
accomplished  after  each  flutter  redesign  cycle  to  account  for  interaction  between 
flutter  and  strength  requirements. 

Results  of  redesign  studies’ using  the  two  cantilever  - structure  models 
and  the  free-free  model  are  presenutd  and  discussed  below. 


12.2  DISCUSSION  OF  REDESIGN  RESULTS 


1 Two  aerodynamic  design  load  conditions  were  apecified  for  the  intermed- 

i iate  - complexity  wing,  namely,  M *■  0.5,  aea  level,  <*  » 38°,  and  M ■ 0.9,  sea 

| level,  a * 9.3°.  Uniform  initial  gages  (0.16  in.)  were  selected  for  all  100 

elements  of  tns  structures  model,  and  a minimum  manufacturing  gage  of  0.02  in. 
was  specified  for  all  cover,  rib  web,  and  spar  web  elements.  The  structure 
was  resized  for  strength  in  five  FSD  cycles  resulting  in  a weight  reduction  of 
the  idealized  structure  from  IU3.0  lb.  to  6l.l  lb.  Adequate  convergence  to  a 
fully  stressed  design  was  indicated  by  the  weight  change  in  the  fifth  redesign 
cycle,  which  was  only  0.1%  of  the  total  final  weight.  The  final  gage  sizes  of 
1 the  fully  stressed  structure  are  presented  in  Figure  12.1. 

I 

This  initial  application  of  the  strength  optimization  program  terminated 
with  the  computation  of  a flexibility  matrix  for  the  dynamics  model.  Dynamics- 
model  mass  data  corresponding  to  the  selected  degrees  of  freedom  were  then  calcu- 
lated by  applying  realistic  nor.optimum  factors  to  the  weight  of  the  fully  stressed 
j wing-box  structure  and  also  by  separately  accounting  for  the  weight  of  the  over- 

| hanging  structure  not  included  in  the  finite  element  model  (see  figure  11.2).  The 

resulting  total  wing  weight  was  187.3  lbs.  compared  with  6l.l  lbs.  for  the  ideal- 
ized primary  structure. 

Flutter  resizing  was  performed  using  all  100  elements  of  the  structures 
model  as  flutter  redesign  variables.  Flutter  analyses  performed  in  the  flutter 
optimization  program  utilized  the  doublet-lattice  unsteady  aerodynamics  routine 
and  the  p-k  solution  procedure.  The  flutter-critical  flight  condition  was  desig- 
nated as  Mach  0.8,  sea  level. 

The  required  flutter  speed  increment  from  the  fully-stressed  design  point 
was  specified  to  be  30%  with  an  acceptable  flutter  speed  band  of  30-31%.  It  was 
decided  to  achieve  the  desired  flutter  speed  in  four  combined  strength/flutter 
redesign  cycles  (user  input).  In  all  subsequent  redesign  cycles  the  mid-band 
flutter  speed  (30.5%)  was  automatically  selected  as  the  target  value.  The  design 
steps,  indicated  in  Figure  12.2,  show  convergence  to  the  design  optimum  in  eight 
cycles,  although  a design  very  close  to  the  optimum  point  was  achieved  after  only 
four  cycles. 
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Figure  12.3  shows  the  distribution  of  the  total  weight  added  to  flutter- 
critical  elements  in  the  final  design.  Also  shown  are  the  final  gages  of  the 
flutter-critical  elements.  It  is  interesting  to  note  that  the  weight  added  to 
the  second  bay  inboard  from  the  tip  accounts  for  approximately  half  of  the  total 
weight  increment. 

The  level  of  uniformity  of  the  final  flutter-velocity  derivatives  of  the 
critical  members  may  be  observed  in  Figure  12.4,  whi^h  also  shows  the  original 
values  of  these  derivatives  for  the  fully  stressed  design.  The  tabular  summary 
of  final  design  data  (Table  12. l)  indicates  that  only  0.6  lbs.  of  structural  weight 
was  eliminated  from  the  strength-critical  structure  in  the  course  of  redesign. 

Thus  the  resizing  required  for  flutter  had  only  a secondary  effect  on  the  struc- 
ture's internal  load  distribution. 

It  should  be  noted  that  strength  resizing  may  be  optionally  by-passed  in 
any  resizing  cycle  when  it  is  apparent  that  strength-flutter  interaction  is  insig- 
nificant. In  this  particular  calculation,  strength  resizing  was  accomplished  in 
each  cycle  simply  to  demonstrate  the  total  program  capability. 

12.2.2  All-Nfovable  Stabilizer 

All-movable  stabilizer  design  load  distributions  were  computed  for  two 
flight  conditions:  M = 0.8,  sea  level,  or  = 16.5°,  and  M = 1.3,  10,000  ft., 
or  * 9-0°.  and  the  preliminary  gages  of  the  structures -model  were  resized  for 
strength  in  five  FSD  cycles.  Since  the  chordwise  and  spanwise  shear  panel 
elements  in  the  stabilizer  structures  model  simulated  the  stiffness  properties 
of  an  aluminum  honeycomb  core,  they  could  not  be  logically  resized  for  either 
strength  or  flutter  requirements.  Consequently  they  were  eliminated  from  the 
redesign  process  by  setting  their  maximum  and  minimum  manufacturing  gages  equal 
to  their  initial  gages.  The  weight  of  the  idealized  structure  was  reduced  by 
15 i (39-5  lbs.)  in  strength  redesign.  The  initial  strength-optimization  program 
analysis  terminated  with  the  computation  of  a flexibility  matrix  for  the  92 
degrees  of  freedom  selected  for  the  dynamics  model  (see  Figure  11.4). 

The  dyncuaics  model  mass  matrix  was  then  calculated  by  applying  non- 
optimum factors  to  the  weights  of  all  the  finite  elements  in  the  structures  model 
ar.d  then  distributing  these  weights  to  the  selected  dynamic  points.  This  hand 
calculation  was  performed  by  a weights  engineer. 
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Figure  12.4  Comparison  of  Flutter-Velocity  Derivatives  (Knots/Lb) 

of  Flutter-Critical  Elements  for  Intermediate  Complexit 
Wing  - 3efore  and  After  Redesign 
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TABLE  12.1  - FINAL  DESIGN  DATA  FOR  THE 

INTERMEDIATE-COMPLEXITY  WING  AND 
ALL-MOVABLE  STABILIZER  DEMONSTRATION 
PROBLEMS 


WEIGHT  IN  POUNDS 


INTERMEDIATE-  ALL-MOVABLE 

COMPLEXITY  WING  STABILIZER 


Weight  of  finite  element  model  - 
preliminary  sizing  of  members 


143.0 


259-7 


- - - after  FSD 

61.1 

220.2 

Total  weight  of  FSD  structure 
including  non-optimum  factors 

187.3 

414.5 

Total  weight  increment  of 
flutter-critical  elements 
from  FSD  (final  design) 


Total  weight  change  of  strength' 
critical  elements  from  FSD 


l4.?6 


-0.60 


Structure  * 3.71 
Mass  balance  = 6.78 

-0.24 


l 


Total  weight  increment  from 
FSD 


14.16 


10.25 


A supersonic  flutter-critical  flight  condition,  Mach  1.6,  30,000  ft,  was 
selected  for  flutter  redesign  with  the  design  goal  being  a 25$  increase  in  flutter 
speed.  Two  separate  studies  were  then  performed,  one  using  both  structural  and 
mass  balance  design  variables  and  the  other  using  only  structural  design  variables. 
In  the  former  case,  3 Its.  of  initial  mass  balance  was  arbitrarily  added  at  each 
of  three  selected  mass  balance  locations  (see  Figure  12.5).  Reasonably  large 
values  of  mass  were  selected  to  avoid  the  possibility,  observed  in  previous 
studies,  of  mass  balance  being  ineffective  for  low  initial  values  even  though  it 
becomes  effective  for  larger  values. 

The  results  of  the  combined  strength/flutter  redesign  study,  including  mass 
balance  variables,  showed  that  the  mass  balance  at  point  B was  increased  beyond 
its  initial  value  whereas  the  masses  at  points  A and  C were  progressively  elimi- 
nated. It  was  obvious,  after  the  third  flutter  redesign  cycle  that  the  mass 
balance  at  points  A and  C would  be  zero  in  the  final  design;  it  was  decided  there- 
fore to  override  the  existing  program  values  and  set  them  equal  to  zero  at  this 
design  point.  This  had  the  effect  of  accelerating  convergence  to  the  design  opti- 
mum point. 

It  is  seen  in  Figure  12.6  that  a near  optimum  design  with  combined. structural 
and  mass-balance  variables  was  achieved  after  only  five  flutter  redesign  cycles. 

The  net  weight  increase  to  achieve  a 25$  increase  in  flutter  speed  was  only  10.7 
lbs  or  2.58$.  A further  three  redesign  cycles  beyond  this  point  eliminated  an 
additional  0.5  lbs.,  a relatively  insignificant  amount.  The  distribution  of  mass 
balance  and  structural  weight  in  the  final  design  is  presented  in  Figure  12.7. 

The  standard  deviation  of  flutter-velocity  derivatives  for  the  nine  flutter  criti- 
cal elements  in  the  final  design  was  0.33  knots/lb  with  a rean  value  of  17.7  knots/ 
lb.  A summary  of  final  design  data  is  presented  in  Table  32.1  wherein  it  is  noted 
that  resizing  for  strength  (a  single  cycle  of  FSD  in  each  combined  resizing  cycle) 
resulted  in  a very  snail  reduction  in  the  weight  of  the  strength-critical  struc- 
ture. 

The  second  all-movable  stabilizer  study,  which  used  only  structural  de- 
sign variables,  indicated  that  the  flutter  effectiveness  of  the  structural  ele- 
ments in  the  tip  region  was  governed  by  their  :sass- balance  effect,  i.e.,  the 
kinetic-energy-density  contribution  was  dominant.  FASTO?  therefore  resized  these 
elements  to  achieve  a mass  balance  effect  similar  to  that  noted  in  the  previous 
study.  Since  the  initial  values  for  the  mass  of  these  minimum  gage  structural 


173 


as 

<H  C 
O V o 
v *H 

5 S Is 
H7l  8 

ropq 


Figure  12.5  Selected  Locations  for  All-Movable-Stabilizer  Initial  Mass  Balance 
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Firwre  12.6  Results  of  All-Movable  Stabilizer  Redesign  Study- 
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elements  was  very  small,  it  became  apparent  that  convergence  to  the  opeimum  com- 
bination of  "mass  balance"  and  structural  stiffness  would  be  slow  (see  Figure 
12.6).  Consequently,  redesign  was  terminated  before  convergence  was  achieved. 

This  case  illustrates  a fact  noted  in  previous  studies  using  FASTOP,  namely  that 
the  convergence  for  flutter  cases  that  are  susceptible  to  mass  balance  is  much 
enhanced  if  the  user  initiates  redesign  with  arbitrary  selections  of  mass  balance 
in  potentially  effective  regions  of  the  structure. 

To  provide  a further  base  of  reference  for  the  two  studies  described  afcovcj 
a final  analysis  was  performed  using  mass  balance  alone  at  the  most  effective 
location  indicated  by  FASTOP,  i.e.,  the  tip  leading-edge.  It  was  determined  that 
14.7  lb.  of  mass  balance  would  be  required  to  achieve  the  25/6  flutter  speed  in- 
crease, compared  with  10.25  lb.  for  the  combined  redesign  case  (see  Figure  12.6). 

12.2.3  Wlng-wlth-Store 

Conversion  of  the  cantilever-structure  dynamics  model  of  the  wing-with- 
store  (see  Figure  11.6)  to  a free-free  model  was  achieved  by  including  beam 
elements  to  simulate  the  stiffness  properties  of  the  fuselage  in  vertical  bending. 
The  fuselage  modelling  is  schematically  illustrated  in  Figure  12.8  in  which  it  is 
noted  that  fifteen  elements  were  used  to  represent  the  fuselage  and  ten  beam  ele- 
ment node  points  were  selected  as  dynamic  node  points  for  vibration  analysis. 

Thus  the  modified  vibration  model  had  a total  of  158  dynamic  degrees  of  freedom 
plus  three  plug  ( rigid  body)  degrees  of  freedom. 

The  preliminary  gages  of  the  finite  element  structures  model  were  resized 
for  strength  in  five  FSD  cycles  using  two  computed  subsonic  aerodynamic  design 
load  distributions  and  one  store  inertial  load  condition.  The  weight  of  the  re- 
sulting fully-stre3sed  structure  was  1340  lbs.,  based  on  the  finite  element  ideal- 
isation, and  1921  lbs,,  including  non-optimum  factors  and  overhangirg  stricture. 

6 2 

The  pylon-mounted  store  weighed  4500  lbs.  with  a pitch  inertia  of  8 x 10  lb.  in 
about  its  center  of  gravity.  The  high  store  inertia  created  a critical 
flutter  mechanism  involving  the  first  wing  bending  mode  and  the  store- pitch/ wing- 
torsion  mode.  A Mach  0.8,  sea  level  flutter-critical  flight  condition  was  de- 
signated for  the  redesign  study.  The  objective  was  to  achieve  a flutter  speed 
target  of  660  knots  equivalent  airspeed  from  an  initial  computed  value  of  270  knots 
for  the  fully  stressed  design.  Wing  posts  and  fhselage  beam  elements  were  excluded 
from  strength/ flutter  redesign  resulting  in  a total  of  453  structural  design 
variables  for  this  study. 
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Figure  12.8.  Ving-with-Store  Dynamic  Model  .of  Fuselage 


T!;e  forward  beam  of  the  pylon  support  structure,  connecting  :.o-i»s  l”f  vA. 
1“7  (See  Figure  11. 5)  had  the  highest  flutter  velocity  derivative  in  the  initial 
design  (236  knots 'lb. ) and  the  beats  element  which  continued  the  forward  rear 
ccnreetion  through  the  lower  and  upper  wing  covers  (nodes  176-175)  had  the 
second  highest  derivative  (97  kr.cts/lh.).  The  rear  beam  of  the  pylon,  which 
had  a flutter  velocity  derivative  of  -hi.  5 knots/lb.,  was  not  resized  downward 
because  of  minimum  gage  constraints.  As  redesign  proceeded,  the  front  and  rear 
beams  of  the  pylon  became  cnually  flutter-effective.  ..'.is  change  was  accorne*- 
ied  by  a noticeaole  change  in  the  flutter  mode  shape,  whlc.i  initially  involved 
components  of  the  first  wing  bending  mode  and  the  store -pitch/wing- torsion  mode, 
but  later  exhibited  increasingly  large  components  of  the  second  wing  bending 
mode  as  the  store  pitch  frequency  was  increased.  A net  reduction  in  structural 
weight  was  achieved  in  the  first  combined  flutter-strength  redesign  cycle 
(Figure  12.9)  due  to  the  fact  that  the  initial  fully  stressed  design  was  not 
fully  converged  (i.e.,  the  stress  ratios  were  not  unity)  and  the  subsequent 
strength  redesign  further  reduced  the  weight  of  the  structure  in  regions  which 
were  ineffective  for  flutter.  Thus,  the  net  weight  reduction  of  O.95  lb.  in 
the  first  redesign  cycle  comprised  an  addition  of  0.3^  lb.  for  flutter  speed 
improvement  and  a reduction  of  1.29  lb.  in  the  strength-critical  regions  of 
the  structure. 

Redesign  for  flutter  appeared  to  be  directed  toward  increasing  the  fre- 
quency of  the  store  pitch  mode  and  no  structural  elements  were  resized  outboard 
of  the  wing  store  station.  Resizing  of  wing  structural  elements  inboard  of  the 
store  station  involved  spar  webs  and  the  rib  webs  between  the  wing-to-pylon 
connection  points  (Figure  12.10).  The  resizing  in  the  vicinity  of  the  front 
wing- to- fuselage  connection  point  accounted  for  a relatively  small  proportion 
of  the  overall  weight  increase.  One  of  the  most  interesting  results  of  this 
study  was  that  the  shear  rigidity  of  tha  forward  and  rear  spars  was  shown  to 
be  a more  significant  contributor  to  the  overall  store  pitch  stiffness  the:! 
the  wing  covers. 

Convergence  to  the  final  design  point  was  achieved  in  eight  combined 
flutter-strength  redesign  cycles  (figure  12. 9),  although  a design  verv  close 
to  optimum  was  achieved  after  six  cycles.  The  flutter  speed  target  was  achieved 
fcr  or.ly  8.^2  lbs  increase  in  weight.  The  distribution  of  this  weight  is  tab- 
ulated below: 
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Front  and  Rear  Wi 
Connection  Points 


Total  weight  increase  for  flutter-criticsl  eleirents  * 
(Weight  increase  for  pylon  alone  * 
Weight  removed  frcm  strength-critical  structure  * 
Total  weight  increase  from  FSD  starting  point  * 


10.35  lb. 
6.55  lb.) 
1.93  lb. 
8.fe2  lb. 


Although  the  net  reduction  in  weight  of  the  strength-critical  structure  was 
relatively  insignificant,  a redistribution  of  cover  material  was  noted  in  the 
inboard  wing  sections.  Consequently  strength  resizing  (l  FSD  cycle)  was 
accomplished  in  all  but  two  of  the  redesign  cycles.  In  the  other  two  cycles, 

SOF  was  simply  used  to  recompute  a revised  flexibility  matrix  for  vibratJon 
mode  calculations. 

A second  redesign  study  was  initiated  using  both  structural  design  variables 
and  three  mass  balance  design  variables  located  at  dynamic  nodes  6,  9>  ari<3  55 
(see  Figure  11.6).  Twenty  pounds  of  mass  balance  was  arbitrarily  added  at  each 
dynamic  node  before  initiating  redesign.  The  flutter  speed  cf  the  fully  stressed 
design  increased  from  270  knots  to  290  knots  as  a result  of  this  mass  addition. 

It  was  noted  that  the  initial  flutter-velocity  derivatives  of  all  masses  were 
very  close  to  zero.  Consequently  they  were  eliminated  in  the  subsequent  re- 
design analysis. 
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